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PREFACE 

j The  following  report  represents  the  first  semiannual  review 

^ of  research  conducted  under  the  auspices  of  the  Associate  Joint 

Services  Electronics  Program  at  the  Institute  for  Electronic  Science 
at  Texas  Tech  University.  Specific  topics  covered  include,  fault 
analysis,  computer-aided  design,  stochastic  control  and  estimation, 

I 

I decentralized  control,  mathematical  system  theory,  optical  noise,  and 

I pattern  recognition. 
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ABSTRACT 

Research  covering  several  aspects  of  the  fault  analysis  problem  for 
electronic  circuits  and  systems  is  reviewed.  This  includes  basic  mathe- 
matical studies  into  the  concepts  of  fault  diagnosis,  fault  predictions 
and  self  testing.  These  ideas  are  then  combined  to  formulate  an  archi- 
tecture for  an  approach  to  built-in  testing. 

Introduction 

The  fault  diagnosis  problem  may  sub-divide  into  three  fundamental 
areas;  fault  detection,  fault  diagnosis,  and  fault  prediction.  These 
areas  are  closely  coupled  and  all  three  are  required  for  any  practical 
test  system.  The  precise  manifistation  of  any  one  being  dependent  on 
the  specific  context  in  which  the  test  system  is  designed  to  operate. 

For  instance,  in  an  off-line  (periodic)  test  system  the  fault  diagnosis 
algorithm  must  be  able  to  cope  with  components  which  have  deviated  a long 
way  from  their  nominal  values  whereas  in  an  on-line  (,built-in)  test  sys- 
tem one  may  assume  that  a failure  is  spotted  as  soon  as  the  components 
begins  to  drift  away  from  nominal.  As  such,  a fault  diagnosis  algorithm 
for  an  off-line  test  system  must  be  designed  to  cope  with  significant 
nonlinearities  whereas  an  algorithm  for  on-line  testing  may  be  formulated 
in  terms  of  a linearized  model. 

In  the  sequel  the  state  of  our  research  in  fault  diagnosis  and  fault 
prediction  is  reviewed.  These  ideas  are  then  used  to  formulate  the  archi 
tecture  for  a built-in  test  (BIT)  system.  Of  course,  the  mathematical 
theory  derived  for  the  fault  diagnosis  and  fault  prediction  problems  is 
also  applicable  to  off-line  test  systems,  test  systems  designed  for  use 
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at  the  end  of  an  assembly  line,  etc.  The  evnironment  in  which  such  test 
systems  operate  is,  however,  quite  different  then  that  of  the  BIT  systems 
and  hence  the  mathematical  techniques  for  fault  diagnosis  and  fault  pre- 
diction would  have  to  be  integrated  differently  in  these  other  applications. 
Fault  Diagnosis 

For  the  purposes  of  doing  fault  diagnosis  we  work  with  a component 
connection  model  for  the  circuit  or  system  under  test  which  takes  for  form 


1. 

b^  = Z.(s,r)a. 

and 

a = L^^b  + L^2 

2. 

y = L2^b  + L22 

6 12  12 

in  the  frequency  domain  ’ ’ . Here  Z^. (s,r)  is  the  transfer  function  of 

the  ith  circuit  or  system  component  where  r = col(r^)  is  the  vector  of 
unknown  component  parameters  and  s is  the  complex  frequency  variable.  The 
L^.j  are  known  connection  matrices,  a = col(a^)  and  b = col(b.j)  are  com- 
posite vectors  of  component  inputs  and  outputs  respectively,  and  u and  y 
are  the  test  input  and  output  signals  respectively.  In  the  nonlinear  case 
the  component  equations  are  replaced  by  the  state  models 

3.  X.  = f^(X.,a.,r) 

i i— 1,2,  ...  ,n 

= g^(X.,a.,r) 

with  the  connection  equations  remaining  as  in  2.  Although  these  component 
connection  models  for  a circuit  or  system  are  nonclassical  they  are  widely 
used  in  large-scale  system  simulation  and  computer-aided  circuit  design  and 
are  readily  amenable  to  the  "computer  speed-up  techniques"  developed  for 
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12 

these  applications.  As  such,  they  are  ideally  suited  for  the  fault 
diagnosis  problem. 


Combining  1.  and  2.  yields  the  fault  diagnosis  equation.^ 

4.  s'"  = L22  + l2i(l-Z(s.r)L^^)"^Z(s.r)L^2 

where  Z(s,r)  = diag(Z^ (s,r) ) and  s'"  is  the  measured  transfer  function 
relating  the  input  test  signal  u to  the  output  test  signal  y.  The  solu- 
tion of  the  fault  diagnosis  problem  therefore  amounts  to  the  solution  of 
* for  the  parameters  vector,  r,  given  s'"  and  the  connection  matrices. 

‘hough  it  is  possible  to  give  an  analytic  description  of  all  possible 

12  13 

solutions  to  this  equation  ’ given  any  fixed  value  for  the  complex 
frequency  variable,  s,  in  a "real  world"  situation  the  number  of  unknowns 
greatly  exceeds  the  number  of  equations  and,  as  such,  the  analytic  repre- 
sentation of  the  solution  manifold  proves  to  be  of  little  value.  This 
difficulty  is  alleaviated  via  a multi -frequency  diagnosis  algorithm  where- 
in one  writes  the  set  of  simultaneous  equations 

5.  S(s.|,r)  = L22  + L2i(l-Z(s^,r)L^l)'^Z(s.|,r)L^2 

S(s2,r)  - L22  i-21  ( i ”Z ( S2 , r ) L-j  1 ) Z(s2,n)Li2 

S(S|^,r)  = L22  ^21  ^ ^ ^ ^1 1 ^ 

where  k different  complex  frequencies  are  used  in  equation  4.  simultaneously. 
The  interesting  and  somewhat  surprising  result  is  that  the  additional  equa- 
tions in  5.  may  be  independent  thus  increasing  the  number  of  equations  with- 
out Increasing  the  number  of  its  unknowns.^  While  the  set  of  simultaneous 
equations  5,  often  have  a unique  solution  no  analytic  solution  technique  is 
known  and  we  must  restort  to  time  consunming  numerical  solution  procedures 
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carried  out  off-line. 

Although  the  multi -frequency  fault  diagnosis  equations  of  5.  do  not 
admit  an  analytic  solution  their  numerical  solution  can  be  significantly 
speeded  up  by  careful  analysis  of  the  equations.  In  particular,  a little 
algebra^’^^  will  reveal  that 


6. 


dS 


m 


= 1 M-Zfs  rlL  '■ /t-11/  j i-ip 

ar  ^iS.,r;L^^  M ii 

^ J 

showing  that  one  can  compute  the  partial  derivaties  required  for  the  nume- 
rical solution  to  5.  analytically.  Moreover,  if  one  observes  that  the 
inverse  matrix  required  to  compute  the  partial  derivaties  in  equation  6. 
is  precisely  the  same  inverse  matrix  required  to  evaluate  the  multi-frequency 
fault  diagnosis  equations  5.  it  is  seen  that  the  partial  derivative  informa- 
tion is  obtained  at  virtually  no  computational  cost  over  and  above  that 
required  for  the  evaluation  of  the  equations.  In  a similar  vain  one  can 
reduce  the  computation  required  to  compute  the  inverses  at  different  complex 
frequencies  by  intergrating  the  differential  equation 


7. 


d(1-Z(s.r)L„)-’  , (i-Z(s.r)L.,r'  hi  C-Hs.DL,,)-’ 

ds  II  as 


using  the  inverse  computed  at  one  particular  frequency  as  a starting  point. 
Although  of  extremely  high  dimension  this  equation  is  easily  integrated  with 
out  the  requirement  for  any  matrix  inversions.  With  the  aid  of  these  obser- 
vations it  is  therefore  possible  to  carry  out  an  entire  interation  of  a 
Newton-Raphson  algorithm  for  the  solution  of  the  multi -frequency  fault  diag- 
nosis equations  with  the  aid  of  only  a single  matrix  inversion. 

Although  one  does  not  have  a "neat"  set  of  equations  such  as  those  des- 
cribed above  for  the  solution  of  the  fault  diagnosis  problem  in  a nonlinear 


2.14 
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circuit  or  system  surprisingly  similar  computational  techniques  can  be 
invoked  in  the  nonlinear  case.  The  key  to  these  techniques  is  the  re- 
placement of  the  mul ti -frequency  information  of  the  linear  case  by  a 
family  of  integral  performance  measures  on  the  test  signals,  u and  y. 

These  play  exactly  the  same  role  in  nonlinear  fault  diagnosis  as  played 
by  the  frequence  information  in  the  linear  case,  allowing  one  to  formulate 

multiple  independent  fault  diagnosis  equations  from  the  same  test  signals. 

3 12 

In  the  nonlinear  case  the  sparse  tableau  algorithm  ’ is  used  to 

evaluate  the  fault  diagnosis  equations  at  each  interation  of  a Newton- 

Raphson  algorithm.  As  in  the  linear  case  this  algorithm  allows  one  to 

compute  the  derivative  information  required  for  the  Newton-Raphson  algorithm 

with  essentially  no  additionally  computational  cost  over  and  above  that 

3 4 12 

required  for  the  evaluation  of  the  equations.  ’ ’ As  such,  by  optimally 
exploiting  the  computational  efficiencies  inherent  in  the  sparse  tableau 
formulation  for  an  electronic  circuit  or  system  it  is  possible  to  obtain 
significant  computational  gains  in  the  solution  of  the  fault  diagnosis 
equations  in  the  nonlinear  case  as  well  as  the  linear  case. 

Fault  Prediction 

Although  for  any  particular  device  one  can  collect  statistical  data 
on  which  to  base  a fault  prediction  algorithm  in  a practical  setting  where- 
in the  same  fault  prediction  algorithm  is  multiplexed  through  the  testing 
of  many  different  SRA's  it  is  necessary  to  use  an  algorithm  which  is  in- 
dependent of  the  specific  properties  of  the  parameter  under  test.  As  such, 

20 

we  employ  a curve  fitting  algorithm.  Although  less  accurate  than  a statis- 

18  19 

tically  based  algorithm  we  have  shown  by  computer  experiment  ’ that  such 


an  algorithm  can  be  employed  as  a satisfactory  fault  predictor.  Such 
algorithms  are  computationally  simple  thus  permitting  a single  central 
microprocessor  to  be  multiplexed  through  tne  testing  of  a large  number 
of  SRA’s.^^ 

Over  Che  past  several  years  we  have  investigated  several  approaches 

to  the  fault  prediction  problem^’^^’^®’^^’^^’^^.  The  first  is  extremely 

18  19  20 

naive  but  has  yielded  surprisingly  effective  results  in  simulation.  ’ ’ 

Basically,  one  collects  data  at  periodic  intervals,  fits  the  data  with  a 
second  order  polynomial,  and  solves  the  quardradic  equation  to  estimate 
the  time  at  which  the  parameter  will  go  out  of  tolerence.  Although  such  an 
approach  might  at  first  appear  to  be  so  naive  as  to  be  inapplicable,  it  has 
yielded  surprisingly  good  results.^®’^^  This  is  due  to  the  fact  that  one 
is  not  really  interested  in  the  accuracy  of  the  failure  time  estimate  but 
only  the  accuracy  of  the  binary  decision  (based  on  this  estimate)  as  to 
whether  or  not  replace  the  SRA.  The  point  is,  that  this  binary  decision  is 
only  made  when  failure  is  expected  in  the  near  future,  a region  of  time  in 
which  a polynomial  extrapolation  is  reasonably  accurate.  On  the  other  hand, 
if  failure  is  estimated  to  take  place  in  the  distant  future  even  though  the 
polynomial  estimate  may  be  in  significant  error  this  will  not  lead  to  an 
erroneous  binary  replacement  decision.  I.e.,  if  failure  is  estimated  to 
take  place  in  3 years  even  if  the  estimate  is  off  by  90«,  the  decision  to 
not  replace  the  SRA  at  this  time  will  still  be  correct. 

A fault  prediction  algorithm  based  on  the  above  described  second  order 
polynomial  extrapolation  has  been  extensively  studied  by  Tung  and  authc>" 

18  19 

and  some  10,000  complete  operations  of  the  algorithm  have  been  simulated.  ’ 
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For  the  most  part  these  simulations  were  carried  out  on  artifical  data 
generated  by  a library  of  special  functions  to  which  a noise  term  was 
added.  These  special  functions  included  some  highly  complex  non-monotonic 
curves.  Additionally,  curves  based  on  the  empirical  drift  formula  for  thin 
film  resistors  were  studies  (R(t)  = At^  where  a lies  between  .3  and  .5). 

In  both  cases  random  noise  with  amplitudes  of  up  to  25%  of  the  tolerence 
interval  was  added  to  the  data.  The  result  of  these  simulations,  which 
we  believe  to  represent  an  environment  which  is  more  extreme  that  the 
"real  world",  was  that  99.5%  of  all  SRA's  were  replaced  before  on-line 
failure  at  a cost  of  about  10%  of  their  lifetime. 

At  the  present  time  a somewhat  more  sophisticated  fault  prediction 

5 

algorithm  is  under  development.  This  is  still  essentially  a curve  fitting 
algorithm  though  one  in  which  a failure  model  well  founded  in  modern 
reliability  theory^  is  employed.  The  basic  idea  for  this  algorithm  is 
as  follows.  The  drifting  SRA  parameter,  r,  is  assumed  to  satisfy  a 
difference  equations 

8.  r(k+l)  = r(k)  = f(k) 

where  the  "component  time"  k represents  the  number  of  shocks  the  SRA  has 
received;  switching  processes,  electrons  boiling  off  a cathode,  etc.  The 
relation  between  component  time,  i.e.  the  number  of  shocks  received,  and 
"real  time"  is  assumed  to  be  a Poisson  distributed  random  variable  in  which 
the  probability  of  the  SRA  receiving  n shocks  in  a time  interval  of  length 
t is 

9.  P^(t)  = (ct)"e'^V  Id_ 
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For  the  fault  prediction  algorithm  it  is  assumed  that  the  value  of 
the  parameter,  r,  is  known  at  a fixed  set  of  points  in  "real  time"; 
r{t.|),  r(t2),  ...  , Using  this  data  we  desire  to  estimate  the 

unknown  failure  dynamics,  f(k),  for  the  SRA  parameter.  This  is  then  used 
in  equation  8.  to  compute  the  number  of  shocks  required  to  cause  failure; 
i.e.  the  smallest  value  of  k fo>'  which  r(k)  is  out  of  tolerence.  Finally, 
this  estimate  is  used  to  compute  the  optimal  "real  time"  at  which  to 
replace  the  SRA  so  as  to  miniminze  the  cost  functional 

10.  J = + c W 

f f w 

Here,  is  the  probability  of  on-line  failure,  W is  the  average  percentage 
of  SRA  lifetime  which  is  wasted  by  replacing  the  SRA  before  its  actual 
failure  and  c^  and  c^  are  weight  factors. 

Note  that  the  implementation  of  the  above  described  Poission  shock 
based  fault  prediction  algorithm  requires  that  we  deal  simultaneously  with 
two  unknown  phenomena;  the  failure  dynamics,  f(k),  and  the  random  relation- 
ship between  "real  time"  and  "component  time"  given  by  the  Poission  dis- 
tribution. Although  the  required  analysis  is  complex  a surprisingly 
tractable  (and  optimal  in  an  appropriate  sence)  fault  prediction  algorithm 
can  be  formulated.  Here,  one  uses  the  properties  of  the  Poission  distri- 
bution to  estimate  the  number  of  shocks  which  the  SRA  has  received  in  the 
time  intervals  [t.,  t._i];  i=l,  2,  ...  ,m  and  combines  this  data  with  a 
generalized  inverse  algorithm  to  estimate  f(k).  Here,  f(k)  is  approximated 
by  a jth  order  polynomial  and  one  must  compute  the  generalized  inverse  of 
an  m by  j matrix.  Fortunately,  the  algorithm  is  ideally  suited  to  a 
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sequential  least  squares  techniques  and  no  matrix  inversions  need  to  be 
carried  out  on-line.  Once  f(k)  has  been  estimated  to  a satisfactory  level 


of  accuracy  (one  keeps  on  increasing  the  order  of  the  approximating 
polynomial  until  the  estimation  error  is  reduced  to  a prescribed  level) 
it  is  used  with  equation  8.  to  compute  the  number  of  shocks,  required  for 
the  parameter  to  go  out  of  tolerence.  Finally,  this  value  is  used  in 
conjunction  with  the  Poission  distribution  to  determine  the  optimal  "real 
time"  at  which  to  replace  the  SRA.  Although  apparently  complex  this  latter 
optimization  can  be  reduced  by  analytic  techniques  to  the  solution  of  a 

5 

single  nonlinear  equation  in  one  variable.  As  such,  the  entire  fault 
prediction  algorithm  miy  be  easily  implemented,  on-line.  Unlike  the  second 
order  curve  fitting  algorithm  the  Poission  shock  algorithm  for  fault  pre- 
diction is  still  under  development  and  its  simulation,  hopefully  on  "real 
world"  data,  is  just  beginning. 

."elf  Testing 

An  interesting  side  effect  of  running  a fault  analysis  system  in  a 
predicitve  mode  is  that  it  opens  up  the  possibility  of  reliable  self  testing. 

I The  key  observation,  here,  is  that  to  do  fault  prediction  in  a digital  de- 

i vice  one  must  test  its  analog  parameters  such  as  rise  time,  power  supply 

i 

I voltages,  clock  speeds,  pulse  widths,  etc.,  since  the  digital  parameters  are 

I either  right  or  wrong  and  have  no  gray  region  within  which  to  extrapolate 

' trends.  Now,  if  one  uses  a microprocessor  to  predict  its  own  failure  by 

extrapolating  the  values  of  its  analog  parameters  as  long  as  the  prediction 

I 

is  made  before  these  parameters  go  out  of  tolerence  the  digital  performance 
of  the  microprocessor  will  still  be  exact  and  hence  one  may  use  the  micro- 
processor as  a reliable  predictor  of  its  own  failure. 


-9- 


4 


The  point  is  that  in  a predictive  mode  the  microprocessor  is  still 
working  at  the  time  it  predicts  its  own  failure  and  hence  may  be  used 
reliable  in  a self  testing  mode.  Of  course,  once  the  analog  parameters 
of  the  microprocessor  have  exceeded  their  tolerence  limits  it  may  no  longer 
be  trusted  as  a digital  signal  processor  and  hence  the  device  cannot  be 
used  to  diagnose  its  own  faults  after  failure. 

Although  the  above  described  self  testing  concept  is  purely  con- 
ceptual and  has  yet  to  be  implemented  nor  even  simulated  it  is  indicative 
of  the  potential  of  fault  prediction  in  a BIT  system.  Indeed,  if  one  can 
reliably  predict  failure  before  it  actually  takes  place  such  "far  out" 
concepts  as  self  repair  move  into  the  realm  of  feasibility  since  at  the 
time  a replacement  decision  is  made  the  device  under  test  is  still  working. 
A BIT  Architecture 

Our  basic  architecture  for  a Built-in  Test  (BIT)  system  is  a two- 
level  hierarchical  structure  illustrated  in  Figure  1.  Intuitively,  the 
overall  system  may  represent  a printed  circuit  board 
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while  the  subsystems  represent  various  shop  replacable  assemblies  (SRA's) 
such  as  integrated  circuits,  power  supplies,  SCR's  vaccumm  tubes,  etc. 
Alternatively,  the  overall  system  may  represent  an  entire  electronics 
system  with  the  SRA's  being  its  consituant  PC  boards.  In  either  case  the 
SRA's  may  be  throw-away  units  or  units  intended  for  off-line  repair  with 
build-in  test  equipment  (BITE)  designed  to  detect  and/or  predict  faults 
in  the  SRA.  For  those  units  intended  for  off-line  repair  the  BITE  may 
also  be  used  as  an  interface  with  an  external  test  stand  but  will  not  be 
capable  of  isolating  the  failure  within  the  SRA. 

This  structure  is  motivated  by  the  above  described  research  into  the 
relative  computational  complexity  of  the  three  fundamental  problems  of 

g 

fault  analysis;  fault  detection,  fault  diagnosis,  and  fault  prediction. 

18  19 

The  latter  problem  requires  considerable  computational  power  ’ but  need 
only  be  carried  out  at  widely  spaced  test  intervals,  say  one  test  per  hour 
(minute,  second,  ?).  As  such,  a single  central  microprocessor  can  be  time 
multiplexed  through  the  testing  of  a large  number  of  subsystem  parameters 
thereby  achieving  the  required  computational  power  for  the  fault  prediction 
algorithm  while  still  holding  the  amount  of  dedicated  test  equipment  within 
reasonable  bounds. 

While  fault  diagnosis  can  be  carried  out  with  considerable  success  the 
process  requires  significant  computational  power  (at  least  a mini  by  today's 
standards)  and  lenghty  computer  runs^’^^’^^.  As  such,  fault  diagnosis  with- 
in an  SRA  is  done  off-line  or,  an  external  test  stand  containing  the  required 
mini  (or  maxi)  computer.  Each  SRA,  however,  will  include  sufficient  BITE, 
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say  a nanoprocessor,  to  collect  and  condition  test  data  on  the  SRA  to  be 
periodically  communicated  to  the  central  microprocessor  for  purposes  of 
fault  prediction  and  detection. 

Fortunately,  both  of  these  endeavors  may  be  achieved  using  a linearized 

model  of  the  SRA  about  its  nominal  values  and  hence  can  be  implemented  with 

12 

relatively  little  computational  power  built  into  the  SRA.  In  particular, 
for  fault  prediction  one  is  interested  in  tracking  various  internal  para- 
meters of  the  SRA  as  they  drift  from  nominal  to  their  tolerence  limit. 

Since  the  tolerence  interval  is  typically  only  a few  percent  this  can  be 
achieved  with  a linearized  model.  For  castastrophic  errors  a linearized 
model  may  be  used  to  detect  failures  even  though  it  is  not  sufficiently 
accurate  for  fault  diagnosis.  As  such,  the  BITE  within  an  SRA  may  be  kept 
within  reasonable  bounds  while  still  delivering  sufficient  data  to  the 
central  microprocessor  for  its  fault  prediction  and  fault  detection  tasks. 

If  needed  fault  diagnosis  within  an  SRA  will,  however,  be  done  off-line 
with  the  BITE  simply  serving  as  in  interface  between  the  SRA  and  an  ex- 
ternal test  stand. 

A final  aspect  of  the  BIT  architecture  is  the  communication  link  be- 
tween the  SRA's  and  the  central  microprocessor.  Here,  one  desires  to  keep 
the  wiring  between  the  SRA's  and  the  central  microprocessor  at  a minimum 
and  simultaneously  would  like  to  have  all  data  transmitted  to  the  central 
microprocessor  in  a uniform  format  so  as  to  permit  interchangability  of 
component  parts  within  the  system.  Although  the  details  of  this  communi- 
cations link  have  yet  to  be  formalized  the  existence  of  an  active  computing 
capability  in  each  SRA  gives  one  considerable  flexibility.  As  such,  we 
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believe  that  is  will  be  possible  to  work  with  a single  test  busJ^ 
the  central  microprocessor  requests  data  from  the  individual  SRA's  by 
transmitting  a signal  on  the  bus.  This  signal  is  received  by  the  built-in 
nanoprocessor  in  the  SRA  which,  in  turn,  transmits  appropriately  condi- 
tioned test  data  back  to  the  central  microprocessor  on  the  same  bus. 

The  above  described  BIT  architecture  would  seem  to  achieve  most  of 
the  requirements  for  a built-in  testing  system. 

i).  Continous  on-line  fault  prediction  and  detection,  up  to  an 
SRA,  is  achieved. 

ii).  The  system  includes  an  interface  for  off-line  fault  diagnosis 
within  an  SRA. 

iii).  Dedicated  test  equipment  represents  a small  percentage  of  the 
total  system. 

iv).  Busing  is  minimized  and  test  data  is  transmitted  to  the  central 
microprocessor  in  a uniform  format  thereby  facilitating  component 
interchangability. 


-13- 


r 


References 


1.  Barlow,  R.E.,  and  F.  Prochan,  Statistical  Theory  of  Reliability  and 

Life  Testing:  Probability  Methods,  New  York,  Holt,  Rinehard  and 

Winston,  1975. 

2.  Chao,  K.-S.,  and  R.  Sacks,  "Integration  Methods  in  Circuit  Analysis", 
IEEE  Proc.,  (to  appear). 

3.  Hachtel,  G.D.,  Brayton,  R.K.,  and  F.G.  Gustavson,  "The  Sparse  Tableau 
Approach  to  Network  Analysis  and  Design",  IEEE  Trans,  on  Circuit  Theory, 
Vol.  CT-18,  pp.  101-113,  (1971). 

4.  Hachtel,  G.D.,  and  R.A.  Rohrer,  "Techniques  for  the  Optimal  Design 
Synthesis  of  Switching  Circuits".  IEEE  Proc.,  Vol.  55,  pp.  1864-1877, 
119671: 

5.  Lu,  K.-S.  Ph.D.  Thesis,  Texas  Tech  University.,  Lubbock,  Texas.,  (in 
preparation) . 

6.  Ranson,  M.N.,  and  R.  Saeks,  "The  Connection  Function  - Theory  and 
Application",  Int.  Jour,  of  Circuit  Theory  and  its  Application,  Vol. 

3,  pp.  5-21,  (1975). 

7.  Ranson,  M.N.,  and  R.  Saeks,  "A  Functional  Approach  to  Fault  Analysis", 

in  Rational  Fault  Analysis.  New  York,  Marcel  Dekker  Inc.,  (to  appear). 

8.  Rao,  C.R.,  and  S.K.  Mitra,  Generalized  Inverse  of  Matrices  and  its 
Appl ications.  New  York,  J.  Wiley  and  Sons,  1971. 

9.  Saeks,  R.,  etal.,  "Fault  Analysis  in  Electronic  Circuits  and  Systems", 

Tech.  Rpt.,  Texas  Tech  Univ.,  Lubbock,  Texas.,  Nov.  1975. 

10.  Saeks,  R. , "An  Experiment  in  Fault  Prediction  II",  Proc.  of  the  1976 

IEEE  AUTOTESTCON,  Arlington,  Tx.,  Nov.  1976,  p.  53,  (abstract  only). 

11.  Saeks,  R.,  Singh,  S.P.,  and  R.-w.  Liu,  "Fault  Isolation  via  Components 
Simulation",  IEEE  Trans,  on  Circuit  Theory,  Vol.  CT-19,  pp.  634-640, 
(1972). 

12.  Saeks,  R.,  and  R.A.  DeCarlo,  "Interconnected  Dynamical  Systems". 

Marcel  Dekker  Inc.,  (to  appear). 

13.  Saeks,  R.,  Wise,  G.,  and  K.-S.  Caho,  "Analysis  and  Design  of  Inter- 
connected Dynamical  Systems",  in  Large-Scale  Dynamical  Systems.  No. 
Hollywood,  Point  Lobos  Press,  1976,  pp.  59-79. 

14.  Saeks,  R. , and  K.-S.  Chao,  "Continuations  Approach  to  Large  Change 
Sensitivity  Analysis",  IEEE  (London)  Jour,  on  Electronic  Circuits  and 
Systems,  Vol.  1,  pp.  11-16,  (1976). 


-14- 


15.  Sangani,  S.,  Unpublished  Notes,  University  of  Detroit,  Detroit, 

Mi.,  1976. 

16.  Sen,  N.,  M.S.  Thesis,  Texas  Tech  University,  Lubbock,  Texas,  1975. 

17.  Tsui,  Y.-L.,  M.S.  Thesis,  Texas  Tech  University,  Lubbock,  Texas,  1975. 

18.  Tung,  L.,  M.S.  Thesis,  Texas  Tech  University,  Lubbock,  Texas,  1975. 

19.  Tung,  L. , and  R.  Saeks,  "An  Experiment  in  Fault  Prediction",  Un- 
published Notes,  Texas  Tech  University,  Lubbock,  Texas,  1976. 

20.  Tung,  L.,  Liberty,  S.R.,  and  R.  Saeks,  "Fault  Prediction  - Towards 
a Mathematical  Theory",  in  Rational  Fault  Analysis,  New  York,  Marcel 
Dekker,  (to  appear). 


T 


-15- 


t 


» 


» 


RESEARCH 

on 

COMPUTER-AIDED  DESIGN 


' K.S.  CHAO 

DEPARTMENT  OF  ELECTRICAL  ENGINEERING 
TEXAS  TECH  UNIVERSITY 


Abstract 


The  use  of  continuation  methods  in  the  computer-aided  analysis  of 
electronic  circuits  is  surveyed.  Such  methods  are  especially  suitable 
when  one  desires  to  compute  the  solutions  to  a family  of  circuit  analysis 
problems  as  a function  of  a continuous  parameter.  Applications  of  the 
concept  to  the  location  of  multiple  solutions  to  nonlinear  equations, 
the  computation  of  input-output  characteristics  for  nonlinear  networks, 
large-change  sensitivity  analysis,  and  the  computation  of  multivariable 
Nyquist  plots  are  reviewed. 

Introduction 

5 9 12  27  33 

Over  the  past  decade  computer-aided  design  > » » • has  become  a 

fundamental  tool  of  the  electronic  circuit  designer.  Indeed,  the  design 
of  integrated  circuits  containing  several  thousands  of  elements  is  now  a 
routine  procedure  in  the  semiconductor  houses,  yet  one  that  could  not  be 
successfully  carried  out  without  computer  aids.  The  difficulty  here  does 
not  lie  wholely  with  the  complexity  of  the  circuits  involved,  though  that 
is  a major  justification  for  the  use  of  computer  aids,  but  also  with  the 
fact  that  the  parasitic  effects  of  the  integration  process  cannot  be 
readily  included  in  a circuit  breadboarded  from  discrete  components.  As 
such,  the  experience  (and  inertia  ) of  the  circuit  designer  has  slowly 
given  way  to  the  onslaught  of  the  computer. 

This  same  period  has  seen  the  development  of  a number  of  sophisticated 
codes  for  the  analysis  of  large  electronic  circuits  such  as  SCEPTRE,  NET, 
CIRCUS,  etc.,  and  an  equal  number  of  interactive  codes  which  permit  the 
circuit  designer  to  interface  with  his  computer  aids  through  some  sort  of 
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graphic  display.  For  the  most  part  these  codes  were  developed  in  in- 
dustry under  the  pressure  of  a specific  application  and  thus,  except  for 
the  use  of  sparse  matrix  techniques’  (required  to  handle  the  large  sys- 
tems of  equations  which  are  encountered),  they  use  routine  numerical  methods. 

In  the  past  several  years  this  situation  has  begun  to  change  the  ad- 
vent of  new  circuit  theoretic  techniques  specifically  designed  for  im- 
plementation in  circuit  analysis  codes.  Most  notable  of  these  have  been 

9 25  26 

the  sparse  tableau  algorithms  of  Hachtel , Brayton,  ’ ’ and  the  conti- 

nuation algorithms  of  Branin,^’^  and  Broyden,^^  Davidenko,^^’^^  Chua  and 
Ushida,^^  for  the  solution  of  the  sets  of  simultaneous  nonlinear  equations 
arising  in  the  analysis  of  nonlinear  circuits.  The  purpose  of  the  present 
report  is  to  survey  these  continuation  algorithms  and  a number  of  related 
continuation  algorithms  applicable  to  the  problems  of  computeraided  analysis 
of  electronic  circuits.  These  latter  include  an  algorithm  for  the  computa- 
tion of  the  input-output  characteristics  and/or  the  AC  analysis  of  nonlinear 
resistive  networks,  an  algorithm  for  the  large-change  sensitivity  analysis 
of  linear  circuits,  and  an  algorithm  for  the  computation  of  multivariable 
Nyquist  plots. 

The  basic  idea  of  all  continuation  methods  is  to  convert  the  solution 
of  a parameterized  family  of  algebraic  problems,  P(r,),  into  the  solution 
of  a differential  equation. 

X (r)  = f[x(r),  r]  1. 

where  x(r)  is  the  solution  of  the  rth  problem,  P(r).  Then,  if  one  can  find 

#This  literature  has  been  reported  in  several  contexts.  For  the  general 
theory  of  sparse  matrices  see  references,  9,  24,  26,  37,  43,  and  49,  for 
application  to  computer-aided  design  see  refs.  2,  6,  20,  22,  25,  27,  30, 

31,  32  and  49;  and  for  applications  to  power  networks  see  ref.  41,  42, 

44  and  45. 
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the  solution,  x(r^),  of  an  initial  problem,  P(rQ),  by  classical  (or 
other)  methods,  the  solutions  to  the  other  problems  can  be  obtained 
by  integrating  equation  (1)  with  x(r^)  as  an  initial  condition. 

Although  it  appears  that  we  have  made  our  problem  harder  by  con- 
verting the  solution  of  an  algebraic  problem  into  the  solution  of  a 
differential  equation  this  is  not  the  case,  and  in  fact,  when  the  in- 
tegration of  equation  (1)  is  reduced  to  a numerical  process  it  often 

proves  to  be  less  cumbersome  than  direct  methods  of  solution.  Indeed, 

★ 


this  is  illustrated  by  classical  quasi-continuation  algorithm  for  the 
solution  of  a set  of  simultaneous  nonlinear  equations.  Here  we  desire 
to  find  an  n-vector  which  satisfies 

f(x)  =0  2. 

where  f is  a continuously  differential  function  mapping  to  R^. 

Rather  than  solving  (2)  directly,  however,  we  integrate  the  differential 

equation^’7,8,12,23 


x(r)  = f[x(r)]  3. 

which  has  a stable  equilibrium  at  the  solutions  of  (2).  As  such,  if  we 
begin  with  any  initial  condition,  x(r^),  sufficiently  close  to  a solution 
the  trajectory,  x(r),  will  converge  to  the  solution  as  r goes  to  infinity. 
Thus,  upon  reducing  the  integration  of  equation  (3)  to  a numerical  algorithm, 
an  algorithm  for  solving  the  original  set  of  nonlinear  equations  is  de- 
fined in  the  process.  In  particular,  if  Euler  integration  is  used  we  obtain 


*This  is  not  precisely  a continuation  algorithm  in  the  sense  defined  above 
since  the  solution  of  the  given  problem  is  the  limiting  value  of  x(r)  and 
we  really  have  no  interest  in  the  intermediary  values.  The  underlying 
philosophy  is,  however,  the  same. 
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the  iterative  scheme 


x(k+1)  = x(k)  - h (|^)‘^  f[x(k)]  4. 

which  will  be  recognized  as  the  classical  Newton-Raphson  algorithm  for 
the  solution  of  the  given  equation  if  the  step  size  h is  set  equal  to 
one. 

Another  classical  example  of  the  solution  of  a family  of  algebraic 
problems  by  continuation  is  the  inversion  of  a family  of  matrices,  M(r), 
via  the  solution  of 

M"^r)  = -M‘^(r)M(r)M'^(r)  5. 

where  M is  computed  by  classical  means.  Indeed,  the  technique  can 

be  used  to  compute  the  inverse  of  a single  matrix,  N,  by  letting 

M(r)  = rN  + (l-r)I  6. 

in  which  case  M(0)  = I and  no  initial  inverse  need  be  computed.  Once 
again,  however,  when  the  integration  of  the  differential  equation  is  re- 
duced to  a numerical  process  a classical  series  expansion  for  the  matrix 
inverse  results. 

The  point  is,  that  even  though  we  have  taken  a roundabout  approach 
to  formulating  our  algorithms,  the  resultant  numerical  processes  are  no 
more  complex  than  those  which  might  have  been  formulated  directly,  Indeed, 
in  the  above  example  our  algorithms  coincides  with  a classical  algorithm 
for  the  solution  of  the  given  problem.  Moreover,  in  a true  continuation 
algorithm  where  one  is  interested  in  the  entire  trajectory,  x(r),  rather 
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than  only  its  final  value  the  resultant  numerical  processes  often  prove 
to  be  more  efficient  than  classical  algorithms. 

Over  the  past  several  years  such  continuation  algorithms  have  re- 
ceived considerable  interest  from  the  numerical  analysis  community. 
Applications  to  polynomial  root  finding,  boundary  value  problems  in 
ordinary  differential  equations,  parameter  identification  and  eigen- 
value problems  in  differential  operators  are  reviewed  in  the  recent 

48 

paper  of  Wasserstrom  and  will  not  be  discussed  here.  Other  significant 
applications  of  the  continuation  concept  include  the  solution  of  Lp 

4 

optimization  problems  by  continuation  with  the  classical  solution 

taken  as  an  initial  condition  and  the  systematic  location  of  multiple 

7 8 13 

solutions  to  sets  of  simultaneous  nonlinear  equations.  ’ ’ 

Multiple  Solutions  of  Nonlinear  Equations 

7 8 

A systematic  search  method  based  on  Branin's  approach  ’ has  recently 

13  14 

been  developed  by  Chao,  Liu  and  Pan  ’ for  obtaining  multiple  solutions 
of  a nonlinear  equation 

f(x)  = 0 7. 

where  f is  a continuously  differentiable  function  from  into  itself. 

The  method  is  capable  of  finding  all  of  the  solutions  provided  that  the 
solution  of  any  (n-1)  equations  of  (7)  is  a simple  curve  - a continuously 
differentiable  curve  which  does  not  intersect  itself.  This  technique 
serves  as  the  key  to  the  development  of  several  of  the  continuation 
algorithms  to  be  discussed  and  is  hence  reviewed  here. 
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The  algorithm  involves  numerical  integration  of  a set  of  associated 


differential  equations  of  the  form 

^^■[x(t)]  = -f^-[x(t)],  f.[x(0)]  = 0,  i = 1,  2,  ...,n-l 

8. 

f,[x(t)]  = .yx(t)],  yx(0)]  = f^^ 

or  in  the  x-space 

X = (3f/3x)'^  (-f,  ....  -f^_^  ^0^  ^ 9. 

along  the  space  curve,  l,  defined  by  the  intersection  of  solution  mani- 
folds for 

f^[x(t)]  =0  i = 1,  2,  ....  n-1  10. 

The  transition  in  sign  of  f^  should  be  made  at  the  solution  points  and 
points  where  the  Jacobian  determinant  changes  sign.  Equation  (9)  may  be 
solved  by  any  numerical  integration  technique.  For  example,  using  Euler's 
method,  (9)  reduces  to  an  iterative  algorithm 

Vl  = \ ^ hJ-’(x^)[-f,(x^), 

x^eA,  k = 0,  1 , 2 11 . 

It  is  interesting  to  note  that  with  only  the  minus  sign  and  for  any  initial 
Xq  not  necessarily  on  t,  (10)  reduces  to  the  well  known  Netwon"s  method. 

It  is  observed  from  (8)  that  if  x^  zi,  the  corresponding  trajectory  x(t) 
resulting  from  (9)  remains  in  l.  Since  in  computational  practice,  the  f , 
for  i = 1,2,...,  n-1,  at  each  iteration  using  (9)  will  not  be  precisely 
zero,  their  signs  are  kept  at  the  constant  negative  to  insure  that  the  com- 
puted trajectory  does  not  stray  too  far  from  i.  If  i is  simple,  we 
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can  obtain  all  the  solutions  by  a complete  traversal  of  the  space  curve 
2,.  In  applying  (9)  to  locate  all  the  solutions,  it  is  thus  essential 
that  there  exists  such  a simple  curve.  In  addition,  a starting  point 
lying  on  or  close  to  this  space  curve  must  be  used.  Theorems  that 
guarantee  this  existance  of  a simple  curve,  t,  and  the  global  conver- 
gence from  any  initial  guess  to  a point  on  i are  all  detailed  in  14  and 
will  not  be  repeated  here.  The  algorithm  is  not  only  capable  of  finding 
multiple  solutions,  it  also  efficiently  computes  the  solution  of  any 
(n-1)  equations  of  (7),  the  distinct  property  of  which  is  found  useful 
in  many  engineering  applications. 

Input-Output  Characteristics 

Operating  points,  driving  point  plots  and  transfer  characteristic 

plots  are  basic  concepts  of  fundamental  importance  in  the  analysis  of 

nonlinear  resistive  networks.  The  operating  point  problem  is  nothing 

but  the  determination  of  the  solutions  of  network  equations  of  the  form 

described  in  (7).  Mathematically  speaking,  the  latter  two  concepts  are 

essentially  the  same;  they  are  input-output  characteristic  plots  showing 

the  relationship  between  a driving  source  u^.  and  a certain  output  vari- 

16  29  34 

able  X..  Several  methods  ’ * are  available  for  piecewise-1 inear 

analysis  of  resistive  nonlinear  networks.  Recently,  Chua  and  Ushida^^ 
developed  a switching-parameter  algorithm  for  solving  nonlinear  equations 

and  the  technique  is  also  capable  of  finding  input-output  plots. 

15 

A technique  is  proposed  for  obtaining  input-output  characteristic 
plots  for  nonlinear  resistive  networks  where  the  characteristic  curve  of 
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each  nonlinear  element  is  continuously  differentiable.  The  network 
equation,  in  general,  can  be  represented  by 

f(x,u)  = 0 12. 


where  f is  a continuously  differentiable  function,  x is  an  n-vector  of 
network  variables,  the  m-vector  u denotes  the  input  sources  and  the 
Jacobian  matrix  3f/3x  is  assumed  to  be  nonsingular  at  the  solution 
points.  Since  in  finding  input-output  characteristic  plots,  all  input 
sources  except  one,  say  u^ , are  considered  to  be  fixed  constrants,  the 
network  equation  (12)  thus  reduces  to 


f(x,u.)  = 0. 


13. 


Suppose  it  is  now  required  to  find  the  input-output  characteristic  plot 
X vs.  u^,  for  all  uT  ^ u^  1 li|.  Instead  of  solving  (13)  directly,  a 
set  of  (n+1)  equations  consisting  of  (13)  and  an  auxiliary  equation 
f^^l  = 0 of  the  form 


f(x,  u.)  = 0 


n+1 


“I  - "i  * 0 


14. 


is  considered.  Following  the  technique  described  above,  a continuous 
input-output  curve  can  be  traced  automatically  by  integrating  a set  of 
associated  differential  equations 


^ = -f 
dt 


f[x(0),  u(0)]  = f(x  , u.) 


= 0 


df 


15. 


dt^  ^ -Vr  Vl^°^  " ‘^i  ■ ‘^i  “ ^n+1,0 


where  x is  a solution  of  (17)  subject  ot  u.j  = u.. 


t 
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The  initial  conditions  on  f and  f^^-j  are  such  that  the  starting 
point  must  lie  on  the  solution  manifold  of  (13)  which  is  the  desired 
TC  plot.  Since  (f,  ) is  not  a function  of  t explicitly,  the  chain 

rule  of  differentiation  is  applied  so  that  in  the  (x,  u^.)  - space,  (15) 
reduces  to 


= J-1 

-f  ~ 

1 

x(0) 

= 

X 

1 

: u. 

L ’ . 

iVi_ 

_u,  (0)_ 

where 

11  ll 

8X  3U^ 

J = ' 17. 

0 1 

Equation  (16)  may  be  solved  by  any  existing  numerical  integration  tech- 
nique and  the  sign  change  should  again  be  made  at  the  Jacobian  singularities. 
From  the  solution  of  (15) 

f[x(t),  u.(t)]  = f[x(0),  u.(0)]e'^  =0 

’ ’ 18. 

ViCui(0]  = = Vl,0®-^  * 

it  is  seen  that  for  any  [x(0),  u^(0)]  et,  the  corresponding  trajectory 
[x(t),  u^(t)]  resulting  from  (16)  remains  in  The  input-output  charac- 
teristic curve  denoted  by  i can  thus  be  traced  from  u^.  to  ut  automatically. 
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The  advantage  of  this  method  is  that  the  input-output  plot  is  not 
required  to  be  either  input  or  output  controlled.  The  only  restriction 
is  that  the  plot  has  to  be  a simple  curve. 

AC-Resistive  Network  Analysis 

Resistive  networks  containing  one  or  more  time-varying  sources  are 
known  as  ac-resistive  networks.  In  this  case,  the  network  equations  can 
always  be  written  in  the  form 

f[x(t),  u(t)]  = 0 19. 


where  u(t)  is  an  m-vector  of  known  time-varying  sources  and  the  n-vector 
X represents  network  variables.  The  operating  point  of  the  network,  which 
is  the  solution  of  (19),  is  now  a function  of  t.  The  presence  of  the  time 
parameter  greatly  increases  the  complexity  of  the  problem.  Theroetical ly, 
the  problem  can  be  treated  by  simply  solving  the  problem  at  different 
instants  of  time  which  would  probably  be  a very  time-consuming  practice. 

A continuation  algorithm  is  proposed  for  finding  the  solution  of  (19)  for 
all  t automatically. 

Consider  a system  of  differential  equations  of  the  form 


f = -f.  f[x(0),  u(0)]  = 0 


20. 


u = w,  u(0)  = u^ 

where  u is  assumed  to  be  continuously  differentiable  and  its  time  derivative, 
denoted  by  w,  is  a known  function  of  t.  The  initial  condition  on  f is  such 
that  x(0)  = x^  must  be  a solution  point  of  (19)  subject  to  u(0)  = u^. 


In 


the  x-space  (20)  is  equivalent  to 


u = -J‘V  w,  x(0)  = Xq  21. 

where  J ” 3f/3x. 

Equation  (21)  may  now  be  solved  by  using  any  numerical  integration 
technique.  If  J is  never  nonsingular  at  the  solution  points  (as  assumed), 
then  lJ|  will  never  vanish  along  the  trajectory  of  (21).  Therefore  once 
an  initial  solution  point  x^  is  found,  the  solution  of  (19)  for  all  t can 
be  obtained  automatically  from  (21). 

If  f[x,  Uq]  = 0 has  multiple  solutions,  then  all  the  initial  solution 
points  must  be  used  as  initial  conditions  for  (21)  in  order  to  obtain  a 
complete  family  of  solution  curves. 

A more  direct  approach  may  be  obtained  by  differentiating  (19)  with 
respect  to  t yiedling 


dt 


i * ft  u . 0 

3x  3u 


22. 


or 


X . w.  23. 

The  only  difference  between  the  two  proposed  schemes  (21)  and  (23) 
is  the  appearance  of  the  correction  term  -J'V  in  (21).  Theoretically 
along  the  solution  curve  of  (24),  f = 0.  In  actual  computation,  however, 
f may  not  be  identically  zero.  The  use  of  the  minus  sign  in  front  of  f 
in  (22)  tends  to  prevent  divergence  and  hence  to  minimize  computational 
error. 
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Large  Change  Sensitivity  Analysis 


The  sensitivity  of  system  performance  is  an  important  aspect  of 
system  design.  The  system  characteristics  may  vary  with  environmental 
changes,  such  as  temperature  or  uncertainties  in  the  parameters  that 
characterize  the  system.  The  analysis  of  a system's  sensitivity  to 
parameter  variations  is  concerned  with  the  determination  of  changes  in 
network  or  system  behavior  which  will  be  produced  by  changes  of  one  of 
its  parameters.  The  study  of  system  sensitivity  to  small  parameter 
variations  around  a reference  position  reduces  to  the  well  known  per- 
turbation problem.  For  large  parameter  variations,  however,  it  is  an 
entirely  different  problem  which  cannot  be  treated  by  perturbation 
methods. Two  approaches  for  large  change  sensitivity  analysis^^ 
in  two  different  contexts  are  discussed  below. 

Sensitivity  to  large  parameter  variations,  or  global  sensitivity, 
for  nonlinear  resistive  networks  can  be  formulated  in  a manner  similar 
to  the  one  of  finding  input-output  characteristic  plots. 

Suppose  it  is  now  desired  to  find  the  operating  point  of  the  net- 
work as  a function  of  a network  parameter  can  be  written  in  the  func- 
tional form 

f(x,  a.)  = 0 24. 

where  f,  as  before,  is  assumed  to  be  continuously  differentiable.  Since 
(24)  is  in  the  same  form  as  (13),  the  operating  point  x as  a function  of 
can  be  obtained  by  integrating  (16)  with  u^  replaced  by  . 
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40 

Our  second  approach  for  sensitivity  to  large  parameter  variations 

35  36  39 

is  based  on  the  component  connection  model  for  a large-scale  system  ’ ’ 

whose  explicit  algebraic  connection  equations  for  a large-scale  system 
yield  analytic  expressions  for  sensitivity  analysis.  It  is  assumed  that 
the  system  is  characterized  by  a mapping  which  relates  the  output  y to  the 
input  u in  the  frequency  domain  by 


y = Su.  25. 

The  component  connection  model  of  (25)  consists  of  writing  a set  of  linear 


algebraic  equations  of 

the  form 

1 

a 

!•  1 
1 

'-12 

1 y 

i_  ^ 

1 

1 

1 

L^i 

l2^ 

1 


26. 


where  the  L..  are  connection  matrices  and  the  component  input  a is  mapped 

* *J 

into  the  component  output  b by  a linear  transformation 

b = Za.  27. 

It  can  be  shown  by  a straight-forward  algebraic  computation  from  (25),  (26) 
and  (27)  that  S is  given  by 

S = L22  + L2^  (I  - ZL^^)"^  ZL^2.  28. 

With  the  assumption  that  the  output  is  not  a function  of  the  input  explicitly, 
i.e.,  L22  = 0»  (28)  reduces  to 

S = 1-2iQL^2 
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where  Q is  an  intermediate  matrix  defined  by 

Q = (I  - 30. 

In  order  to  derive  an  algorithm  for  computing  the  transfer  function 
matrix  sensitivity  to  large  component  variations,  we  differentiate  (30) 
with  respect  to  a potentially  variable  parameter  r to  yield 

= Q(r)Z-''(r)^^  [I  + L^^Q(r)].  31. 

Since  for  a fixed  system  connection,  Li.j,  Z(r)  and  hence  Z(r)  are  known 
analytically,  the  right-hand-side  of  (31)  can  easily  be  computed  for  a 
given  r.  We  first  compute  Q(r^)  for  a nominal  value  r^  of  the  parameter 
r by  classical  analysis  techniques  and  then  integrate  equation  (36)  for 
all  r^  r^.  Once  Q(r)  is  known,  then  S(r)  can  be  calculated  directly 
from  (34).  Note  that  if  the  variable  parameter  r is  chosen  to  be  the 
angula>^  frequency  u,  then  the  frequency  response  S(ju))  as  a function  of 
u can  be  computed  in  a similar  manner.  Starting  with  an  initial  condition 
Q(i^q).  it  is  seen  from  (31)  that  the  integration  of  Q(r)  involves  only 
matrix  multiplications  and  a simple  inversion  of  the  usually  elementary 
component  matrix  Z(r).  Computationally,  the  technique  thus  requires  much 
less  time  than  would  be  needed  to  reanalyze  the  system  which  would  normally 
require  an  inversion  at  each  step  of  a much  more  complicated  matrix  (I-ZL.|.|) 
as  indicated  in  (30). 

Multivariable  Nyquist  Plots 

Although  the  classical  Nyquist  criterion  had  long  been  thought  to  be 
inapplicable  to  multivariable  systems  Barman  and  Katznelson^  have  recently 
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formulated  a Nyquist-like  test  for  such  systems.  Here,  the  classical 
Nyquist  plot  is  replaced  by  a plot  made  up  of  the  eigenvalue  loci  of  the 
(open  loop  gain)  transfer  function  matrix  G(jw).  More  precisely,  for 
each  fixed  value  of  u,  G(jw)  has  n (not  necessarily  distinct)  eigenvalues 
X^(j(jj).  As  such,  if  one  lets  u evlove  one  can  compute  n eigenvalues  loci 
for  G(jw)  by  plotting  the  n functions  x^(u)  as  functions  of  w.  In  general 
the  x^.(o))  are  not  unique  but  they  can  always  be  formulated  as  piecewise 
analytic  functions.^  We  form  a "multivariable  Nyquist  plot"  from  the 
X^(oo)  by  concatenating  the  analytic  segments  obtained  by  plotting  the 
various  X^(u)  to  form  a closed  curve  in  the  complex  plane.  Using  this  plot 
it  has  been  shown^  that  the  multivariable  feedback  system  with  a stable 
open  loop  gain  G(joj)  will  be  stable  if  and  only  if  its  multivariable  Ny- 
quist plot  does  not  encircle  -1. 

In  order  to  apply  the  multivariable  Nyquist  criterion,  it  is  thus 

necessary  to  compute  the  eigenvalue  loci  of  G(u)  as  a function  of  frequency. 

For  a given  frequency,  the  eigenvalues  can  be  calculated  by  using  classical 

techniques.  Since  the  eigenvalues  are  functions  of  frequency,  normally  one 

would  have  to  repeat  the  entire  computational  procedure  for  each  frequency. 

In  the  actual  stability  analysis,  this  repetition  is,  however,  impractical. 

38 

Our  first  approach  to  the  formulation  of  a continuation  algorithm  for 

. 1 

computing  the  eigenvalue  loci  follows  directly  from  that  described  by  Faddeev 
and  Faddeeva^^  and  Van  Ness  et  al^^.  The  eigenvalues  X^.(a))  of  G(juj)  and 
their  complex  conjugates  X^(u))  satisfy 

G( jw)X^ ((d)  = X^((d)X^((d)  i = 1,  2,  ...,  n 32. 
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and 


G (ju)V.((j)  = i - 1,  2,  n 


33. 


where  X^(oi)  and  V^(u)  are  the  corresponding  eigenvectors  of  x^(aj)  and 

★ 

x^(w)  respectively,  and  G (joj)  is  the  complex  conjugate  transpose  matrix 

of  G{jaj).  We  differentiate  (32)  with  respect  to  oj  to  yield 

Ar  dX.  dx.  '^^i 

d;:^•  ^ ^ " dlT  ^i  ^i  dT- 

Taking  the  scalar  product  on  both  sides  of  (34)  with  Vj  gives 

dX.  dx.  dX. 

< ^ X,,  V.  ^ + < G;:-*-,  V,>  = < X.,  V.  > + x.<  V >.  35. 

dtij  i j duj  j do)  1 J i duj  j 


Using  (33),  the  identity 


dX.  dX. 

> = X.< 

Qu*  J J utLi 


and  letting  j = i , we  obtain 


dG  X.,  V.> 
d X , ^ J — 1 1 

^ dw 

du<  <X.,  V.  > 


, i “ 1,  2,  ...,  n. 


In  deriving  a differential  equation  for  X^ , we  let 
dX.  n 

= r d.  . X.,  i = 1 , 2 n.  38. 

aw  1J  J 

In  view  of  the  fact  that  the  eigenvectors  are  unique,  one  may  assume  with- 
out loss  of  generality  that  = 0.  The  coefficients  for  i ^ j are 
obtained  by  forming  the  scalar  product  with  Vj  on  both  sides  of  (37)  and 
using  the  orthogonality  conditions 


< X.,  Vj  > = 0 i j 


to  yield 


ij 


dX. 

- ST 


Substituting  (36)  and  (39)  into  (35),  leads  to 


40. 


V . = -I- 

^ dii.  ’ ^j'  X.  - \ . ' du  ''i 


dG  V 1/ 

< -X.,  Vj  >, 


Combining  (40)  and  (41)  yields  the  desired  expression  for  a.. 
dG  , „ 

»ij  ■ u,  '’‘J- 

In  a similar  manner  we  obtain 


41 


42. 


i'i 

d(,j 


n 

Z 6-dV.,  i = 1,  2, 
j=l  ^ 


43. 


where 


dVi 

^ dlT’  V 


®ii  " ®ij  " <v,,  x,> 


^ j- 


44. 
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Abstract 


Once  a state  model  for  a stochastic  linear  plant  is  obtained,  the 
analytical  aspects  of  designing  a feedback  controller  can  be  conceptually 
dichotomized.  The  first  part  includes;  selection  of  a performance  measure 
that  reflects  a-priori  design  specifications,  selection  of  performance 
indices,  and  selection  of  a controller  via  optimization  of  these  indices. 

The  second  part  includes  selection  of  post  design  performance  measures  and 
the  attainment  of  statistical  or  probabilistic  descriptions  of  these  per- 
formance measures. 

Introduction 

We  will  collectively  refer  to  the  first  part  as  design  and  the  second 
part  as  performance  analysis.  For  certain  classes  of  stochastic  linear 
control  systems  the  performance  analysis  problem  has  beer,  solved;  see  [7] 
and  [14].  The  purpose  of  this  paper  is  to  present  the  formulation  of  a 
complete  set  of  statistics  of  the  integral  quadratic  design  performance 
measure  normally  employed  in  stochastic  linear  control  system  design.  These 
statistical  formulations  are  expressed  explicitly  in  terms  of  dynamical 
variables  related  to  the  plant  observations  and  control  action.  The  ultimate 
value  of  these  formulations  is  that  they  provide  the  analytical  basis  for 
the  selection  of  design  indices  which  can  be  used  to  systematically  select 
a controller  via  optimization. 

There  are  many  ways  of  selecting  indices  and  we  do  not  treat  this  topic 
here.  In  Section  II  we  describe  the  system,  the  performance  measure  and  the 
control  objective.  Section  III  contains  the  devf'lopment  of  a complete  str..is- 
tical  description  of  this  performance  and  Section  IV  contains  the  formulation 
of  these  statistics  in  terms  of  the  filtered  estimate  of  the  plant  state. 

The  paper  is  concluded  with  comments  on  how  the  results  might  be  employed 
and  suggestions  for  further  research. 
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The  System  Description  and  Performance 


Let  denote  the  n-fold  Cartesian  product  of  the  real  line,  and  let 
I denote  the  real  line  interval  [t^,  t^].  We  wish  to  control  the  noisy 
linear  system  described  on  I by 

= A(t)x(t)  + B(t)u(t)  + C(t),  (1) 

and 

z(t)  = C(t)x(t)  + 0(t),  (2) 

where  the  state  x(t)  cR^,  the  control  action  u(t)  cR^,  and  the  observation 
z(t)  cR"^.  The  initial  condition  for  (1),  x(t  ),  is  assumed  to  be  Gaussian 
with  mean 

Xo  = E{x(tg)}  (3) 

and  covariance 

Zo  = E{[x(t^)-x^][x^t^)-xJ]}  (4) 

where  (^)  denotes  matrix  transposition.  The  state  process  noise,  ^(t) 
and  the  observation  noise  are  zero-mean  Guassian-whi te  with 


E{^(t)0^(T)} 

= 0, 

(5) 

E{[x(t^)-xJ^/(t)} 

= 0, 

(6) 

E{[x(t^)-xj0^(t)} 

= 0, 

(7) 

E{C(t)J(T)} 

= H(MMt-T), 

(8) 

E{0(t)0^(T)} 

= Q(t)6(t-T), 

(9) 
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where  5(t)  and  0(t)  are  positive  semi -def ini te  and  positive  definite,  respectively, 
on  I. 

We  require  that  the  control  action,  u(t),  be  a causal  function  of  the  ob- 
servation. That  is, 

u(t)  = i/^(t,  z(T):Tc[tQ,t])  , (10) 

where  tj/(t,  •)  satisfies  certain  technical  assumptions  stated  in  [14];  also 
see  [15].  All  matrix  functions  on  I and  the  mapping  4i(t,«)  are  assumed  to 
be  smooth  enough  to  guarantee  mean-square  continuity  of  the  state  process  on  I. 

For  the  purposes  of  design  we  define  a "design-performance-measure", 

>1,  by 

J = x^(t.)Sx(t,)  + / ^ [x^(t)Q(t)x(t)  + u^(t)R(t)u(t)]dt,  (11) 

^0 

where  the  terminal  penalty  weighting,  S,  is  symmetric  and  positive  semi- 
definite  as  is  the  weighting  Q(t)  on  I.  The  weighting  R(t)  is  symmetric 
and  positive  definite  on  I,  and  both  Q(t)  and  R(t)  are  continuous  on  I. 

These  weighting  matrices  are  given  values  by  the  designer  that  reflect 
a-priori  design  specifications  involving  the  relative  importance  of  state 
regulation  and  control  effort.  The  design  objective  is  to  choose  u(t) 
in  (10)  so  system  performance  is  "good"  in  some  sense. 

The  functional,  J,  assigns  a non-negative  real  number  to  each  sample  run 
of  the  control  system  with  small  values  implying  good  performance.  However, 
the  question  of  quality  of  performance  is  multiply  clouded.  First,  J is  random 
so  it  is  only  meaningful  to  refer  to  J in  a statistical  or  probabilistic  sense. 
Second,  since  J is  the  sum  of  terms  representing  measures  of  state  regulation 
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and  control  effort,  the  individual  quality  of  these  measures  is  not  apparent 
in  a broad  statistical  description  of  J.  We  will  demonstrate  these  subtleties 
in  another  paper.  Meanwhile  we  are  content  to  use  J as  an  a-priori  performance 
measure  upon  which  we  will  base  our  selection  of  control  action,  u(t).  But, 
before  this  selection  can  be  made  we  must  describe  J statistically  or  prob- 
abilistically. Based  upon  intuition  gained  in  our  previous  work  [7]  we  choose 
a statistical  approach. 


A Complete  Statistical  Description  of  J 


Let  be  the  sigma-algebra  induced  by  the  observation  {z(x) :Tc[t^,a] }. 
When  a = t^  we  shall  write  F without  a subscript.  We  .vill  now  generate  con- 
ditional statistics  of  J . Expand  the  process  modeled  by  (1)  in  an  orthonormal 
series. 


x(t)  ~ E X (})  (t), 
i=l  ^ ^ 


where  the  x^'s  are  scalar  random  variables  given  by 


X.  = x''’(t^)  Sc}>.(t^)  + / x'^(t)Q(t)4..(t)dt,  Vi, 

^0 


and  the  nonrandom  vector  valued  functions,  4>^(t),  are  chosen  to  satisfy  the 


orthonormality  condition 


4>{(t^)S4»j(t^)  + ;%^(t)Q(t)4)j(t)dt  =6.j  , Vi,  j, 


In  addition,  we  require  that  the  x^'s  be  conditionally  uncorrelated,  that  is 


E{[x^-m^][Xj-mj]|F}  = 

where  m^  is  the  conditional  mean  of  x^  given  by 
m^  = E{x^(t^)|F}  S4>^(t^) 


Vi,  j. 


+ / E{x^(t)lF}Q(t)<^.(t)dt, 

^0 


Vi,  j. 


A necessary  and  sufficient  condition  for  (15)  given  (13)  and  (14)  is 
that  4>^(t)  and  satisfy 


^^4>.(t)  = / r(t,x)Q(T)4)  (i)dT 

t.  ^ 


r(t.t^) 


S4>^(t^),  tcl,  Vi, 


(17) 


where  r(t,x)  is  the  smoothed  estimate  error  covariance  kernel  of  the  state 
process.  That  is,  let  the  smoothed  estimate  of  x(t)  be  denoted  by 


^(tlt^)  = E{x(t)|F),  tcl. 


(18) 


and  let  r(t,x)  be  given  by 


r(t,x)  = E{[x(t)  - x(t(t^)]  [x(x)  - x(x(t^)]lF)  t.xcl.  (19) 


As  a consequence  of  the  1 inear-Gaussian  assumptions  of  Section  II  and  the 
technical  assumptions  on  <//(t,*). 

r(t,x)  = E{r(t,x)}  (20) 


Therefore  each  X^  is  nonrandom. 

Under  the  assumptions  we  have  made,  J is  finite  with  probability  one; 
see  Doob  [1].  It  follows  that  the  series  in  (12)  converges  in  the  square 
integrable  sense.  That  is. 


’•f 

J = ? x^  + / u^(t)R(t)u(t)dt, 
i=l  i t 


(21) 


where  convergence  is  with  robability  one;  see  Kolmogorov  and  Fomin  [5  ]. 
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Since  x(t)  is  conditionally  Gaussian  each  is  conditionally  Gaussian. 

But,  we  have  forced  the  x^  to  be  conditionally  uncorrelated,  thus  they  are 
conditionally  independent  as  are  their  squares.  The  conditional  characteristic 
function  of  each  x?  term  in  (21)  is  of  the  noncentral  chi-square  type  given  by. 


C^2|p(jw)  = (l-2j«j^)''^exp[jwm^^(l-2jw^)‘^]. 

The  conditional  characteristic  function  of  J follows  as, 

C,|j-(jw)  = [ S (l-2wX  )"^]  • exp[ju?  u^(t)R(t)u(t)dt 

i=l  i t^ 

+ ? jiuTi^(l -2jtcX^)"^]. 


(22) 


(23) 


f 


In  our  previous  work  [6  ],  [ 7 ] we  have  observed  that  in  the  Linear-Quadratic- 
Gaussian  class  of  systems  the  second  characteristic  function  generates  the 
most  tractable  statistical  forms.  The  second  conditional  character! Stic 
function,  Tjjj-(jw),  is  defined  as  the  natural  logarithm  of  Cjjp(jw),  that  is, 


Tj|p(jw)  = ln[Cj|p(jio)]. 


(24) 


The  MacLaurin  series  representation  of  Tj|p(jw)  is 

i! 


(25) 


where  the  coefficients  K^jpare  called  conditional  cumulants.  Utilizing  (23), 
it  can  be  easily  shown  that  the  first  conditional  cumulant  of  J is  given  by 
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(26) 


'llF  " Jl  ih  u^(t)R(t)u(t)dt, 

^ ^0 


while  the  remaining  conditional  cumulants  are  of  the  form 


k-1 


^klF 


+ k!2'^‘^  mh  , k > 1 

i=l  i 1=1  1 i 


(27) 


Although  the  conditional  cumulants  as  given  by  (26)  and  (27)  are  complete  in 
the  sense  that  any  statistic  of  J can  be  derived  from  them  (an  interesting 

exercise)  they  are  not  in  an  attractive  form  for  the  control  system  designer 

since  they  are  not  expressed  in  terms  of  system  variables.  To  accomplish  this 

we  must  attack  the  series  expressions  in  (26)  and  (27).  The  first  step  is  to 

define  "iterated  kernels" 


and 


(t.x)  ^ r(t.T) 


(t.T)  ^ r(t,  t^)  sr^'^'^^t^.x) 


(28) 


+ ^ r(t,o)  Q(o)r^*^"^  ^ (o,x)  do,  k > 1. 


(29) 


It  can  be  inductively  shown  using  (14)  and  (17)  that 


r^'^^t.x)  = x|j(^^(t)(|j(x). 


(30) 


It  follows  that  the  expression,  ili  written  as 

xj  = Tr[Sr^‘^^t^,t^)  + / Q(t)r^'^^t,t)  dt],  k^l. 


(31) 
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where  Tr  [•]  denotes  the  trace  of  the  enclosed  matrix.  Utilizing  (16),  (18) 
and  (30)  it  follows  that 


= i^(t^|t^)  (t^.t^)  sJ(t^|t^) 

+ i^(t^|t^)S  ^ r^'^'^)(t^,t)Q(t)^(t|t^)  dt 
^0 

+ / ^(tlt^)Q(t)  r^'^"^^(t,  t^)dt  Sx(t^lt^) 

^0 

+ V J^(t|  t,)Q(t)r^'''^)(t,T)Q(T)^(TitJdTdt,  k>i.  (32) 
t t ^ ^ 


For  the  case,  k=l  , it  is  easily  seen  that 

= x^(t^|t^)Sx(t^|t^) 

+ / x (tlt,)Q(t)x(t|t,)dt.  (33) 

t ^ ^ 

0 

Consider  the  last  term  in  (32)  and  note  that  it  contains  a symmetric  (in 
argument)  integrand.  Therefore  it  may  be  rewritten  as 
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(34) 


V-T 


t t ^ ^ 


0 0 


rr(k-T)/ 


= 2/x  (t|t,)Q(t)  /r^  ‘ Mt,T)Q(T)x(T|t,)dTdt 

t ^ t ^ 

0 0 


Define  the  new  n-vector  valued  variable 


_-j(t)  - / ^ ^(t,T)Q(T)x(Tlt^)dx 


and  note  that 


= 0- 


The  conditional  cuxiulants  can  now  be  written  as 


‘1|F 


X^(t^lt^)  Sx(t^lt^)  + /[X^(tlt^)Q{t)X(t!t^) 

^0 


+ u (♦■  )R(t)u(t)]dt 


+ Tr[sr(t.,t,)  + / Q(t)r(t,t)dt], 
^0 


and 


t 
I 
t 


.(k-1). 


+ 2j^t^lt^)Sn,^.l(t^)  + 2 / x^(t|t^)Q(t)n,^_-,(t)dt] 

^0 

+ (k-l)!2'''''Tr[Sr^''^t^.,t^)  + Q(t)r^''^t,t)dt] , k > 1 


(35) 


(36) 


(37) 


(38) 

« / 
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Equations  (37)  and  (38)  provide  us  with  the  conditional  cumulants  of  the 
design-performance  measure,  J,  expressed  in  terms  of  the  smoothed  estimate 
of  the  process  and  the  corresponding  error  covariance  kernel.  A complete 
statistical  description  of  J is  now  at  hand  since  (as  mentioned  ■■'arlier) 
any  statistic  of  J can  be  expressed  in  terms  of  the  conditional  cumulants. 
For  example,  denoting  the  unconditional  cumulants  of  J by  have 


E(J}  = 

(39) 

Var(J}  = <2 

= E{'^2|F^  * Var{K^  ,p} 

(40) 

and  in  general 

tCk  = E{<k|p}+  {statistics  of  lower  order  conditional  cumulants}  (41) 

The  relationship  between  noncentral  moments,  cumulants  is  well-known 

[ 4 ] and  given  by 

^k+1  " j^Q^^^k-j'^j+l  . 
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The  Filtered  Estimate  Formulation 


( 


I 


t 

1 

I 

\ 


Those  familiar  with  the  traditional  minimum  mean  LQG  problem  may 
be  a little  suspicious  of  equation  (37)  since  it  is  well  known  that 
E{tCj^|p)  is  normally  expressed  in  terms  of  the  filtered  estimate  and 
its  corresponding  error  covariance  with  precisely  the  same  structure 
as  (37)  under  expectation;  see  [14].  To  demonstrate  the  equivalence 
of  these  two  formulations  note  that  in  view  of  (10)  the  smoothed 
estimate  can  be  expressed  as 

^ /V  T -1 

x(t|t,)  = x(t|t)  + r K(t,r)  C (t)9  (t)  C(t)  v(Tjr)dT,  (43) 

^ t 

see  [ 8], where 

= E{x(t)|F^},  (44) 

K(t,x)  = Ef[x(t)  - x(t|t)][x^(x)  - x^(x|x)]}  (45) 

and  the  "innovation"  [2],  v(*|*)»  is  given  by 

v(t]t)  = C(t)[x(t)  - x{t\t)-]  + e(t).  (46) 

The  smoothed  error  covariance  can  also  be  expressed  as 

r(t,x)  =K(t,x)-/^  K(t,o)C^(o)'3 (o)  C(o)|((a,x)do  (47) 

tvx 

where  tvx  means  max  [t,x]. 


-50- 


Substitution  of  (43)  and  (47)  into  (37)  and  application  of  expectation 
immediately  yields 


= E{j'^(t^|t^)S^(t^|t^)  + }^[J^(tlt)Q(t)^(t|t) 


+ u^(t)R(t)u(t)]dt}  + Tr[SE(t.)  + / Q(t)Z(t)dt] , (48) 

^ t 

0 


where 


E(t)  5 K(t,t) 


(49) 


and  we  have  utilized  the  fact  thatv(*I*)  is  white  noise  with  covariance 
E{v(t|t)v^(xlT)}  =0(t)6(t-x).  (50) 

In  deriving  (48)  from  (37)  only  one  subtlety  arises  that  might  be 
troublesome  to  the  reader.  In  particular,  two  terms  of  the  form 


E{  x^(tlt)Q(t)  K(t,x)C^(x)0"^ (x)C(x)v(xlx)dt) 

^0  ^ 

= Tr[  / Q(t)  /K(t,x)C^(x)0'^  (x)C(x)E{v(x|x)x''’(tlt)}dxdt]  (51) 

^0  ' 

arise. 

It  is  well  known  [14]  that  under  some  technical  assumptions  on  the  causal 
mapping  \p(t,  •)  in  (10) 

t t 

x(t]t)  = ; G(t,x)v(x  |x)dx  + / <l>(t,x)B(x)u(x)dx  + 4>(t,t  )x  (52) 

0 0 

for  some  G(t,-)«  where  <J>(t,*)  is  the  transition  matrix  associated  with 
A(t)  in  (1).  Since  the  control  u(-)  is  assumed  to  be  a causal  function 
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of  the  observation  z[')  which  in  turn  can  be  expressed  as  a causal 
function  of  the  innovation  v(  • |*  )j  [ 2 ] , it  foil  ows  in  view  of  (50) 
that  the  inner  most  integrand  in  (51)  is  zero  almost  everywhere. 
Consequently  the  terms  in  question  vanish  under  expectation. 

For  higher  order  cumulants  things  aren't  as  nice.  Define  the 
following: 

t 

n (t)  ^ V^.T)Q(T)x(T|i)dT,  (53) 

k-1  ^0 

P (t)  ^ K(t,t)C^(T)  0‘^(T)v(xlT)dT, 

t 

and 

P^;_l(t)^.r  r^'"‘^^(t,T)Q(T)p  (i)dT.  (55) 

^0 

Then  substitution  of  (43)  into  (35)  and  (38)  yields 

= k!2''-''[  x^(t^lt^)  sr^'''''^(t^,t^)s.<^(t^|t^) 

+ 2^^t^|t^)Sfi^_^(t^)  + 2:?^(t^lt^)Sp^_.,(t^) 

+ 2 ^(t|t)Q(t)n,^_^(t)dt  + "^(tlt)Q(t)p^_^(t)dt 

to  tQ 

+ 2 /^p  ^t)Q(t)fi^_^(t)dt  + 2/^  p ^(t)  Q(t)p^_^(t)dt] 
^0  ^0 

+ (k-1)!  2'^'^  Tr[sr^^^t.,tf)+  / Q(t)r^'^^(t,t)dt],k>l . 

^ t 

‘'0 

(56) 
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The  variable  l(t),  defined  in  (53)  can  be  generated  by  a physically 
realizable  linear  dynamical  system  with  the  filtered  estimate  as  input.  Con- 
sequently it  is  implicitly  affected  by  the  control  action  u(t).  This  system 
is  of  order  2n  and  the  derivation  of  its  dynamical  equations  is  left  as  an 
exercise.  On  the  other  hand,  the  variable,  p(t),  defined  in  (54)  cannot  be 
produced  by  a causal  system  but  by  a noncausal  linear  dynamical  system  of 
order  n operating  on  the  innovation  process.  Thus,  p(t)  is  unaffected  by 
the  control  action,  u(t),  as  is  Pj^  i(t)  in  (55). 

Consequently  the  conditional  cumulants  as  given  by  (56)  contain  four 
distinct  types  of  terms: 

i.  terms  independent  of  the  dynamical  variables  and  thus 
uncontrollable, 

ii.  terms  containing  only  causally  produced  variables  that 
are  controllable, 

iii.  terms  containing  only  noncausally  produced  variables  that 
are  uncontrollable,  and 

iv.  terms  containing  both  causally  produced  variables  and 
noncausally  produced  variables. 

The  conditional  cumulant  descriptions  that  we  have  derived  form  an  ana- 
lytical basis  for  the  selection  of  performance  indices  which,  in  turn,  can 
be  utilized  in  the  optimal  selection  of  a controller.  This  dynamical  vari- 
able formulation  is  attractive  for  such  optimization  but  carries  with  it  a 
"curse  of  dimensional i ty"  in  that  system  order  will  increase  linearly  with 
the  number  of  higher  order  pieces  of  statistical  information  included  in 
any  selection  of  performance  indices.  The  noncausally  produced  variables 
are  troublesome.  At  this  point  it  is  not  clear  how  they  should  be  handled 
in  optimization  while  enforcing  (10).  Such  questions  are  the  subjects  of 


future  research. 


V.  CONCLUSIONS 


How  might  these  formulations  be  utilized?  There  are  many  possible 
answers  to  this  question.  For  example,  let  us  assume  that  the  design 
objective  is  to  select  a controller  that  will  keep  the  variance  of  J small 
without  the  mean  of  J becoming  too  large.  Such  an  objective  suggests  that 
we  select  a weighted  sum  of  the  indices,  mean  and  variance,  for  optimiza- 
tion. Thus  we  might  choose  as  our  criterion 

m^n  + a(E{K2|  pJ'’’Var{<j  I p})j 

subject  to  the  obvious  dynamical  constraints.  The  nonnegative  real  para- 
meter, a,  here  is  simply  a design-tradeoff  parameter  between  mean  and  vari- 
ance. Clearly  the  expressions  included  in  Var{Kj^|p}  are  rather  complicated 
(fourth  order)  and  it  might  be  worthwhile  to  drop  this  term  with  the  hope 
that  minimizing  the  parameterized  tradeoff  between  E{ic^jp}  and  Ei<'2|F^ 
will  achieve  the  design  objective  for  some  value  of  a. 

Future  directions  of  research  are  clear.  We  should  attack  the  ques- 
tions of  index  selection  and  optimization.  These  questions  are  fairly 
difficult,  particularly  in  the  continuous  time  format.  Perhaps  the  dis- 
crete time  formulation  which  has  not  yet  been  derived  will  be  more  trans- 
parent. 

We  feel  that  the  major  observation  that  we  have  made  is  that  cumulants 
(in  this  case  conditional  cumulants)  retain  the  second  order  nature  of  the 
original  performance,  J,  whereas  other  statistics  do  not.  The  structural 
simplicity  of  the  conditional  cumulants  is  the  key  to  our  success  in  the 
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derivation.  The  simple  fact  that  other  investigators  have  not  considered 
cumulants  explains  why  the  formulations  have  gone  undiscovered  despite 
many  years  of  vigorous  research  activity  in  the  area  of  Linear  Quadratic 
Gaussian  control  theory. 
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Abstract 


A study  directed  towards  the  development  of  a special  pur- 
pose microprocessor  architecture  for  application  in  decentralized 
control  is  reviewed.  Specific  emphasis  is  placed  on  assembly 
language  routines  for  implementing  vector-matrix  operations  on 
the  microprocessor. 

Introduction 

.Many  equations  of  modern  control,  digital  filtering,  Kal.man 
filtering  and  related  problems  involve  matrix-vector  multipli- 
cations which  must  be  calculated  to  obtain  their  solution.  These 
matrix-vector  products  have  the  general  form  shov;n  in  Equation  1 


A X = 


^11 

^12 

• • 

Si 

Im 

^21 

^22 

• • • 

2m 

^nl 

n2 

a 

nm 

•'ll 

X-  1 

2 

• 

X 

m i 
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where  A is  a n x m matrix  and  x is  a m x 1 column  vector.  In 
computing  a solution  to  problems  where  the  product  A x is  needed, 
considerable  time  is  involved  for  any  reasonable  size  n and  m. 
Hence  if  a real-time,  on-line  solution  to  the  problems  of  modern 
control  and  signal  processing  is  to  be  computed,  a method  for 
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calculating  A x efficiently  must  be  developed.  The  method 
developed  is  to  design  an  appropriate  architecture  and  firm- 
ware which  will  yield  an  efficient  implementation  of  m^atrix- 
vector  product  algorithms.  The  procedure  used  to  derive  this 
appropriate  architecture  and  firmware  is  outli.ned  below. 

First  a machine  language  implementation  of  the  algorithm 
to  compute  A x is  produced  so  that  comparisons  of  execution 
times  with  subsequently  developed  microcoded  versions  can  be 
made.  An  algorithm  is  developed  to  calculate  execution  times 
for  any  n,  m,  and  T where  n and  m refer  to  the  matrix  size  and 
T is  the  basic  m.achine  language  instruction  execution.  Next, 
various  size  matrix-vector  products  are  microcoded.  A n x 2 
matrix  times  a 2 x 1 vector  microprogram  is  implemented  and 
run.  A Hewlett  Packard  HP  2100,  minicomputer  was  used  to 
determine  both  machine  language  and  microcode  execution  times. 
Algorithms  are  developed  to  calculate  execution  times  and  actual 
execution  ti.mes  for  various  size  matrix- vector  products  are 
computed  for  comparison.  Finally  half-word  length  (8-bit)  ver- 
sions of  a n X 2 matrix  times  a 2 x 1 vector  and  a 4 x 4 m.atrix 
times  a 4 X 1 vector  are  microcoded  and  run.  Algorithms  for 
execution  time  calculation  are  developed  and  execution  times 
computed  for  various  size  products  to  compare  with  the  afore- 
mentioned machine  language  and  full-word  length  versions.  This 
development  of  micoroutines  reveals  the  architectural  charac- 
teristics necessary  for  an  effective  implementation  of  matrix- 
vector  products  on  microprogram  controlled  computers. 
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Assembly  Language  Routine 


The  assembly  language  routine  for  an  n x m matrix  times  an 
m X 1 vector  is  shown  flowcharted  in  Figure  1.  The  routine 
was  programmed  to  handle  only  integers  that  had  been  scaled 
three  octal  places,  thus  actually  allowing  three  octal  place 
accuracy.  No  checking  for  overflow  or  underflow  was  made. 

The  multiplication  of  integers  is  certiJnly  justifiable  since 
the  main  concern  is  not  to  just  demonstrate  feasibility  or 
viability  of  this  approach  but  to  develop  routines  amenable  to 
real-time,  on-line  systems  applications.  In  such  applications, 
information  into  the  computer  derives  from  an  analog  to  digital 
converter  which  produces  integerized  digital  values  of  the 
analog  signal.  Scaling  is  normally  employed  between  the  com- 
puter and  the  system  both  on  imput  and  output.  Checking  for  an 
overflow  condition  would  also  normally  be  carried  out,  but  the 
action  taken  upon  an  overflow  detection  is  system  dependent, 
e.g.,  an  abort  might  be  necessary,  further  scaling  might  be 
sufficient,  etc.  Therefore,  the  routine  developed  is  suffi- 
ciently general  and  adequate. 

The  elements  of  A,  the  a^^'s  were  stored  sequentially  by 
rows,  i.e.,  first,  then  through  a then  ^2\' 

The  algorithm  used  for  accessing  each  array  element  is  given  in 
Equation  2. 

= (I  - 1)N  -t-  J + Address  ( a^j^)  - 1 


Address (a ^ j) 


(2) 


The  time  of  execution  for  this  routine  can  be  calculated  from 
the  following  equation, 

t^  = 29.5mnT  + lO.SnT  + T (3) 

where  n and  m refer  to  the  matrix  size  and  T is  the  basic 
machine  instruction  execction  time,  ( 1.96  usee  for  the  HP2100). 
Equation  3,  along  with  all  other  equations  in  thl'S.._paper  re- 
lating execution  times,  was  determined  by  actually  summing  the 
times  to  execute  each  instruction  and  noting  the  dependency 
on  m and  n in  processing  the  algorithm.  The  times  for  a 2 .x  2 
matrix  times  a 2 x 1 vector,  a 4 x 4 matrix  times  a 4 x 1 vector, 
and  an  3 X 8 matrix  times  an  8 x 1 vector  which  were  calculated 
for  latter  comparisons  are,  respectively, 

t^2  = 274.40  us 

t . , = 1009.40  us 

44 

^88  ~ 3867.08  us 

.Microcoded  Full-Word  Matri.x- 
Vector  Multiplication  Routines 

The  flowchart  for  a microcoded  multiplication  of  an  n x 2 
matrix  times  a 2 x 1 vector  is  shown  in  Figure  2.  To  save  memory 
references  and  unnecessary  programming,  x^  and  x^  are  first  read 
into  weS  and  stored  in  registers  F and  Q,  respectively.  All 
scratch-pad  registers  are  used:  S2  and  S4  in  the  multiply  sub- 

routine to  temporarily  hold  the  multiplicand  and  multiplier  as 
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are  A and  3 to  hold  the  results;  S3  contains  the  running  sura, 
and;  SI  to  count  n.  Each  a is  retrieved  from  memory  as 

-L  kJ 

needed;  and  the  result  a^^  is  then  stored  in  memory. 

The  Flag  flip-flop  is  used  to  determine  when  this  result  is 
completed. 

The  algorithm  for  calculation  of  time  of  execution  was  de- 
termined to  be 

t^  = 13T  + llOnT  (5) 

where  n is  the  number  of  rows  an  the  matrix  and  T is  the  micro- 
instruction execution  time,  (196  ns  for  the  HP2100) . The  time 
to  execute  this  microprogram  for  a 2 x 2 matri.x  times  a 2 x 1 
vector  was  calculated  from  Equation  5 to  be  45.688  ys.  Com- 
paring this  with  the  time  to  perform  the  same  operation  in 
assembly  language,  a savings  factor  of  six  (6)  or  about  230  ys 
is  accomplished. 

As  was  demonstrated,  all  available  registers  plus  the  Flag 
flip-flop  were  used  to  implement  this  microprogra.m.  Since  only 
microinstructions  and  no  data  (other  then  eight  bit  constants 
stored  in  the  least  significant  eight  bits  of  certai.n  micro- 
instructions) can  be  stored  in  WCS , a severe  limitation,  this 
is  the  maximum  size  matrix-vector  product  that  can  be  impleme.nted 
on  the  HP-2100.  However,  microprograms  of  a 4 x 4 matrix  times 
a 4 X 1 vector,  an  8 x 8 matrix  times  an  8 x 1 vector,  and  an 
n .X  m matri.x  times  an  m x 1 vector  were  flowcharted  and  actually 
microcoded,  but  not  implemented.  These  flowcharts  are  shown  in 
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Figures  3,  4,  and  5,  respectively.  These  routines  require  one 
storage  register  for  each  element  of  the  vector  (4,  8,  and  m) , 
two  scratch  registers  plus  the  A-  and  B-registers  for  the 
multiply  subroutine,  one  scratch  register  for  the  counter  for 
the  square  matrices  and  four  counters  for  the  n x m matrix, 
and  one  register  for  the  running  sum;  or,  8,  12,  and  m + 7 
scratch  registers,  respectively.  As  seen,  these  are  more  than 
are  available  on  the  2100,  a limitation  that  could  be  overcome 
if  data  could  be  stored  and  retrieved  directly  from  WCS. 

The  time  of  execution  algorithms  for  the  first  two  of  these 
microprograms  is  given  as 

2 n-2 

t = 5T  + 7nT  + 59n  T + nT  > (n  - d) 

® d=2 

n-1 

+ 2nT  ^(n  - d)  , (6) 

d=2 

where  n refers  to  the  (square)  matrix  size,  T is  the  micro- 
instruction execution  time,  and  d is  a counter  which  depends 
on  the  number  of  decisions  to  be  made  (n-dependent) . The  times 
of  e.xecution  for  these  programs  are 


(4x4)  t = 200.116  us 
e 

(8x8)  t = 860.244  us 
e 


(7) 


Comparing  these  times  with  those  corresponding  ones  found  using 
the  assembly  language  routine  (i.e.,  200.116  us  vs.  974.12  us 
and  860.244  us  vs.  3867.08  us),  savings  ratio  of  about  4.8  and 
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4.5  are  effected. 


The  algorithm  for  calculating  e.xecution 


time  for  the  n x m matri.x  times  the  m .x  1 vector  micro- 
program IS 

m-2 

t = 2T  + 5mT  + 58n.mT  + 4nT  + nT  ^ (m  - d) 
® d=l 


m- 1 

+ 2nT  ^ (m  - d)  , (8) 

d=2 


Proposed  Architecture 

From  the  preceeding  discussion  concerning  implementation 
of  microroutines  for  matrix- vector  products,  the  limitations 
of  the  HP-2100  for  efficiently  microprogramming  such  problems 
are  evident.  A proposed  architecture  to  alleviate  these  limi- 
tatio.ns  and  allow  an  efficient  implementation  of  an  n x m 
matrix  times  an  m x 1 vector  is  shown  in  Figure  9 and  des- 
cribed in  the  following  paragraphs. 

.^ssuming  that  the  bus  structure  and  microinstruction  format 
remain  fixed,  the  most  severe  limitation  of  the  HP-2100  is  a 
shortage  of  scratch-pad  registers.  By  providi.ng  more  scratch- 
pad registers,  larger  size  matrix-vector  products  can  be  im- 
plemented, greater  computational  versatility  results,  and  the 
execution  speed  of  many  programs  can  be  increased  considerably. 
To  implement  the  n x m matrix  times  m .x  1 vector  product,  m + 7 
registers  are  needed.  These  m + 7 registers  are  shown  in  Figure 
6,  where  m scratch-pad  registers  are  shown  as  "S-bus"  registers 
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the  other  seven  are  shown  as  "R-bus"  registers.  Four  of  the 
"R-bus"  registers  are  labeled  as  in  the  HP-2100,  i.e.,  A,  B, 

Q,  F.  Providing  more  registers  on  the  R-bus  also  results  in 
greater  versatility,  a reduction  in  number  of  microinstructions 
to  implement  a given  function,  and  hence,  a reduction  in 
execution  time.  All  of  these  registers  are  also  assumed  to  be 
general  purpose  registers  and  not  latches  as  the  scratch-pad 
registers  in  the  HP-2100  are  likely  to  develop  a "’•ace"  condi- 
tion if  loaded  while  being  interrogated.  This  limitation  pre- 
vents the  microprogrammer  from  specifying  the  same  scratch-pad 
register  in  both  the  S-  and  T-bus  (in  the  same  microinstruction) 
which  .leads  to  more  microinstructions  than  necessary  if  these 
registers  are  made  general  purpose. 

Figure  7 also  shows  several  additional  five  bit  counters  on 
the  -S-bus.  These  counters  are  not  necessary  to  implement  this 
matrix-vector  product  since  the  m + 7 registers  include  the 
registers  necessary  for  counting.  The  counters  are  shown  to 
indicate  that  several  of  the  m + 7 registers  may  be  replaced 
by  shorter  length  (5  bit)  registers  to  be  used  as  counters - 

Another  limitation  of  the  HP-2100  is  that  only  microin- 
structions can  be  stored  in  and  executed  from  WCH,  that  is  no 
data  can  be  directly  accessed  in  WCS.  The  only  data  available 
in  WCS  is  stored  as  eight  bit  constants  in  the  least  signi- 
ficant eight  bits  of  microinstructions  containing  a "CR"  or 
"CL"  micro-order  in  the  S-bus  filed.  The  "CR"  and  "CL"  micro- 
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orders  direct  the  computer  to  read  the  eight  bit  constants 
stored  in  bits  0-7  of  the  microinstruction  onto  the  least 


(CR)  or  most  (CL)  significant  bits  of  the  register  specified 
in  the  T-bus  field.  The  feature  to  read  and  store  data 
directly  into  WCS  locations  by  microprograms  resident  in  WCS 
would  greatly  enhance  the  capabilities  of  the  machine.  How- 
ever, the  need  for  this  feature  is  mitigated  by  the  addition 
of  sufficient  scratch-pad  registers  but  if  m is  large,  the 
cost  of  such  registers  might  become  prohibitive.  So  there  is 
a trade-off  here,  either  incorporating  as  many  registers  as 
necessary  to  implement  a problem  or  incorporating  several 
additional  registers  and  adding  the  capability  to  read  the 
store  data  into  V?CS.  These  features  would  significantly  in- 
crese  speed  of  execution  of  many  microprograms  and  greatly 
enhance  the  computer's  overall  capabilities — not  only  for 
matrix-vector  products  but  for  a wide  class  of  problems. 

While  the  research  was  implemented  on  a HP- 2 100  minicomputer, 
the  proposed  architecture  could  be  fabricated  as  a single 
microprocessor  clip  using  LSI  techniques. 
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Fig.  6.  Time  of  e:cecution  vs.  matrix  size. 
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Abstrac 


This  research  is  concerned  with  the  generalization  of  fre- 
quency domain  techniques  in  network  and  system  theory  to  non- 
linear, time-variable,  distributed,  and  digital  netv/orks  and 
systems  via  the  techniques  of  Non-Self-Adjoint  spectral  theory. 
Specific  projects  include  a study  of  the  underlying  mathematical 
foundations  of  the  spectral  factorization,  applications  to' gen- 
eralized Wiener  filters  and  optimal  feedback  controllers,  ana 
the  formalization  of  the  concept  of  an  Orlicz  resol tuion  space. 

Introduction 

The  goal  of  our  on-going  research  project  on  "Mathematical 
System  Theory"  is  to  extend  the  frequency  domain  techniques  of 
classical  network  and  system  theory  to  nonlinear,  time-variable, 
distributed,  and  digital  networks  and  systems  via  the  methods  of 
non-self-adjoint  spectral  theory  applied  to  operators  defined  on 
an  abstract  resolution  space.  As  such,  the  research  has  been 
two-fold  in  nature.  One  aspect  of  the  work  deals  with  the  form- 
ulation of  the  operator  theoretic  results  in  resolution  space 
and/or  one  of  classical  L2  spaces  while  the  second  aspect  of  the 
research  is  directed  toward  the  reformulation  of  classical  fre- 
quency domain  results  to  aid  in  their  generalization  and  the  appli- 
cation of  the  resultant  generalizations  in  the  development  of  new 
frequency  domain  techniques. 

This  dichotomas  approach  is  illustrated  by  our  recent  work  on 
the  formulation  of  a generalized  Nyquist  stability  theory.  The 
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first  step  was  the  development  of  a new  proof  of  the  classical 

theorem  which  made  the  homotopic  nature  of  this  theorem  explicit. 

This,  in  turn,  led  to  the  formulation  of  a Nyquist-like  theorem 

for  abstract  nonlinear  operators  on  a resolution  space  wherein  a 

system  is  shown  to  be  stable  if  its  open  loop  gain  is  homotopic 

1 9 

to  the  identity  operator  in  an  appropriate  sense.  this  condi- 
tion has  been  shown  to  coincide  with  the  classical  Nyquist  cri- 

g 

terion  in  the  linear  time-invariant  case  and  is  believed  to  be 
"tight"  in  the  general  case.  Indeed,  most  of  the  classical  suf- 
ficient conditions  for  stability  including  a generalization  of 

g 

the  circle  criterion  can  be  derived  from  the  homotopy  theorem. 

Finally,  the  intuition  derived  from  this  general  theorem  has 

led  to  the  formulation  of  a new,  and  highly  surprising,  Myqu’st- 

like  criterion  for  multi-dimensional  digital  filters  characterized 

6 7 

by  functions  in  several  comples  variables.  ’ These  results  form 
a totality  which,  we  believe,  illustrates  the  essence  of  our 
approach  for  even  though  the  classical  Myquist  criterion  could 

g 

be  derived  v;itnout  an  explicit  call  to  homotopy  theory  the  homo- 
topic formulation  was  the  direct  predecessor  to  the  abstract 
theorem  in  resolution  space  and  the  multivariable  theorem,  both  of 
which  use  homotopic  ideas  in  a highly  non-trivial  manner. 

A second  and  still  incomplete  aspect  of  the  research  has 
been  the  generalization  of  the  spectral  factorization  theory, 
classically  employed  in  frequency  domain  design  techniques,  to 
abstract  operators  defined  on  resolution  space.  We  have  previously 
demonstrated  the  existence  of  a miniphase  factorization  and  ex- 
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hiDited  its  fundamental  relationship  to  the  theory  of  repro- 

20 

ducing  kernel  spaces.  During  the  past  year  a number  of  vari- 
ations on  this  theory  have  been  studied  and,  in  particular,  it 
has  been  shown  that  much  of  the  theory  can  be  extended  to  Banach 

space  in  such  a way  that  the  reproducing  kernel  space  remains  a 
24 

Hilbert  space.  As  such,  we  have  been  able  to  show  that  the 

reproducing  kernel  space  for  a Banach  space  valued  random  variable 

is  a Hilbert  space  and  that  the  scattering  representation  of  a 

network  with  voltage  and  current  vectors  defined  in  Banach  space 

24 

lies  in  a Hilbert  space.  Needless  to  say  all  of  these  results 

have  been  motivated  by  classical  frequency  domain  theory  and  they 

have  led  us  into  the  study  of  the  spectral  factorization  of  rational 

functions  in  several  complex  variables.  This  latter  study,  which 

is  just  beginning,  appears  to  have  significant  protential  both 

theoretically  and  practically.  We  have  already  shown  that  scalar 

rational  functions  in  several  comples  variables  need  not  have  a 

1 4 

rational  spectral  factorization  and  we  are  in  the  process  of 

developing  a design  criterion  for  mul ti -deminsional  digital  filters 

which  will  assure  that  the  resultant  rational  approximation  to 

the  given  insertion  loss  specification  admits  a factorization 

(and  hence  a miniphase  realization).^^ 

Two  other  areas  of  research  which  have  been  undertaken  during 

the  past  year  include  the  extension  of  the  resolution  space  con- 

2 

cept  to  a relativistic  time  structure  and  the  formulation  of  a 
non-standard  model  for  resolution  space.  In  the  former  case  a 
theory  has  been  formulated  for  the  case  of  two-dimensional  special 
relativity  whereas  in  the  latter  case  only  partial  results  have  been 
obtained. 
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Nyquist  Theory:  To  formulate  the  frequency  domain  Myquist  theory  in 


the  classical  and  multivariable  cases  we  require  the  following  notation. 

1.  p"  = •,(z^.Z2’---’^n^  " I’il  - ^ ^ 

is  termed  the  n-variable  polydi sc  in  the  space  of  complex  n-truples  and 

2.  t'’  = z^)  < C";  ;z.j  = 1; 

is  termed  its  distinguished  boundary.  Of  course,  in  C^,  and 

reduce  to  the  usual  1-variable  disc  and  its  boundary.  For  a function,  f 

in  one  complex  variable,  which  is  analytic  in  its  Nyquist  plot  is 

the  image  of  f restricted  to  T^*.  With  this  notation  the  classical  Nyquist 

g 

theorem  has  the  following  statement. 

Theorem:  For  a function,  f,  in  one  complex  variable  defined  as 

above  the  following  are  equivalent. 

i.  f has  no  zeros  in  P^. 

ii.  f has  no  zeros  in  T^  and  INDgf  = 0. 

iii.  f has  no  zeros  in  T^  and  is  homotopical ly  trivial  when 

viewed  as  a mapping  taking  its  values  in  -0. 

Here,  the  index  of  the  complex  valued  function,  f,  is  defined  in  the 
13 

usual  manner  and  the  equivalence  of  ii.  and  iii.  is  the  classical 

13 

result  attributed  to  Hopf  . Although  the  proof  of  the  theorem  can 

g 

be  obtained  as  a direct  corrollary  to  the  classical  argument  principal 
in  which  the  hcmotopic  nature  of  the  result  is  then  implicit  in  the 
proof  of  the  argument  principal  we  believe  that  the  homotopic  \'iew- 
point  is  fundamental  to  the  nature  of  the  theorem  and  thus  have  formu- 


*Actually  it  suffices  for  f to  be  analytic  in  the  interior  of  P^  and 
continuous  on  T . 
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lated  a new  explicitally  homotopic  proof  of  the  theorem.  In  this  proof 

the  necessity  follows  immediately  via  a simple  continuity  argument  while 

the  sufficiency  proof  requires  that  one  carefully  modify  the  domain  of  f 

to  make  it  a covering  map  whence  the  result  follows  from  classical  homo- 
1 3 

type  theory. 

The  generalization  of  the  Nyquist  theorem  to  functions,  f,  analytic 
in  p'^  takes  the  following  somewhat  surprising  form. 

THEOREM:  For  a function,  f,  in  several  complex  variables  defined  as 

above  the  following  are  equivalent. 

i . f has  no  zeros  in  p'^ . 

ii.  f has  no  zeros  in  T^  and  = 0. 

iii.  f has  no  zeros  in  T^  and  f_  is  homotopical  ly  trivial  when 
viewed  as  a mapping  taking  its  values  in  C^-0. 

Here,  f is  the  function  of  one  complex  variable  constructed  from  f 
via  f(z)  = f(z,z,  ...  , z). 

To  our  knowledge  the  theorem  is  sharper  than  any  known  stability  test  for 
multivariable  digital  filters  in  that  the  required  test  is  n-dimensional 
whereas  existing  stability  tests  are  at  least  (n+1 )-dimensional . Indeed, 
in  applications  one  often  knows  a-priori  whether  or  not  f has  zeros  on 
T^  in  which  case  the  remaining  part  of  the  test  is  1 -dimensional . 

The  proof  of  the  multivariable  Nyquist  criterion  will  appear  in 
reference  ^ and  will  not  be  repeated  here.  In  essence,  one  first  decomposes 
the  multivariable  polydisc  into  a continuum  of  single  variate  polydiscs 
and  applies  the  classical  Nyquist  criterion  to  each  disc.  This  then  results 
in  a test  composed  of  a continuum  of  classical  Nyquist  tests  similar  to  that 
proposed  in  reference  6.  Fortunately,  each  of  the  required  plots  can  be 


shown  to  be  homotopic  to  one  of  n Nyquist  plots  and  hence  this  continuum 
of  plots  reduces  to  a finite  number  of  Nyquist  plots.  Finally  the  analy- 
ticity  of  f is  used  to  reduce  the  remaining  n plots  to  the  single  Nyquist 
plot  for  f. 

Cur  final  theorem  on  Nyquist  theory  extends  the  homotopic  ideas 

described  above  to  finite  gain  nonlinear  operators  defined  on  a resolution 
21 

space.  Here,  we  let  Tq  and  be  unbiased,  finite  gain,  causal  operators 
on  a resolution  space  and  we  say  that  Tq  is  homotopic  to  T,  if  there  exists 
an  operator  valued  function,  T,  defined  on  [0,1]  such  that:  T(t)  is  causal 

and  has  a finite  gain  inverse  for  each  t in  [0,1]  T^(0)  = Tq,  and  T(l)  = T-j . 
THEOREM:  Let  T be  an  unbiased,  finite  gain,  causal  operator  defined 

on  a resolution  sapce  which  is  homotopic  to  the  identity  operator. 

Then  T’^  is  causal . 

This  theorem  is  a natural  generalization  of  our  earlier  result  wherein  T”^ 

is  shown  to  be  causal  if  0 is  in  the  unbounded  component  of  the  spectrum 
22 

of  T.  Indecj,  in  that  case  the  required  homotopy  may  be  constructed  via 

2 

a spectral  mapping.  The  proof  of  the  new  theorem  follows  from  a compact- 

28 

ness  (of  [0,1])  argument  quite  similar  to  that  used  in  the  earlier  theorem. 
Although  the  above  theorem  would  seem  to  be  quite  abstract  most  known 

g 

sufficient  conditions  for  stability  follow  readily  from  the  theorem.  For 

1 9 

instance  one  may  derive  the  circle  criterion  by  constructing  a homotopy 
in  two  stages.  First  the  nonlinearity  is  contracted  to  a scalar  with  the 
diameter  of  the  circle  being  such  as  to  assure  the  operator  is  causal  and 
invertable  as  every  state  of  the  contraction.  Then,  in  the  second  stage  of 
the  homotopy  the  remaining  operator  is  deformed  to  the  idenity  using  the 
"non-enci rclement"  hypothesis. 
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Factorization : The  factorization  problem  in  Hilbert  space  is  essentially 


the  problem  of  factoring  a positive  hermitian  operator,  Q,  as 
3.  Q = KK*  (Q  = L*L) 

where  K (L)  has  an  appropriate  causality  structure.  This  is  an  abstracti- 
fication  of  several  problems  commonly  encountered  in  linear  system  theory 
such  as  the  factorization  of  matrix  values  functions  of  a complex  variable, 
the  solution  of  Wiener-Hopf  equations  and  the  solution  of  matrix  Riccatti 
equations.  As  such,  a thorough  understanding  of  this  problem  is  essential 
to  the  extension  of  the  ideas  of  linear  system  theory  to  new  settings. 

In  reference  20  we  proved  the  existence  of  a unique  (up  to  a resolution 
space  equivalence)  miniphase  factorization  for  an  aribtrary  positive  hermitian 
operator  on  a resolution  space  and,  furthermore,  showed  that  the  resultant 
factor,  K,  could  be  represented  as  the  injection  operator  mapping  the  re- 
producing kernel  resolution  space  for  Q into  the  space  on  which  Q is  defined. 

The  main  thrust  of  our  work  during  the  past  year  has  been  the  extension  of 

94 

this  result  to  operators  mapping  a Banach  resolution  space  to  its  dual. 
THEOREM:  Let  (B,P)  be  a (reflexive)  Banach  resolution  space  and  Q 
be  a positive  self-adjoint  operator  mapping  B to  B*.  Then  there 
exists  a unique  (up  to  a resolution  space  equivalence)  miniphase 

operator,  K,  Mapping  a Hilbert  resolution  space,  (H,E),  to  (B,P)  such 

★ 

that  Q=KK  . Furthermore,  K may  be  represented  as  the  injection 
operator  mapping  the  reproducing  kernel  resolution  space  for  Q into 
(B,P). 

24 

The  proof  of  the  theorem  follows  lines  quite  similar  to  the  proof  of  the 
Hilbert  space  theorem  of  reference  20  as  does  the  definition  of  the  re- 
producting kernel  (Hilbert)  resolution  space  for  a positive  hermitian  operator 
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mapping  B to  B. 

The  important  and  surprising  result  of  the  theorem  is  that  even 
when  one  is  working  with  operators  on  a Banach  resolution  space  the 
factor  space  is  assured  to  be  a Hilbert  resolution  space.  This,  in 
turn  implies  that  when  one  applies  such  factorization  theorem  to  systems 
defined  on  a Banach  resolution  space  the  resultant  constructions  will  take 
the  form  of  a Hilbert  space.  Two  such  applications  which  we  have  been  in- 
vestigating are  the  reproducing  kernel  space  wherein  one  may  represent  a 
Banach  resolution  space  valued  random  variable  as  white  noise  (random 
variable  with  identity  covariance)  in  its  associated  reproducing  kernel 
Hilbert  resolution  space.  Similarly,  we  have  shown  that  the  scattering 
representation  for  an  electric  network  defined  in  Banach  space  is  a mapping 

between  two  (possibly  different)  reproducing  Kernel  Hilbert  resolution 
24 

spaces. 

As  indicated  earlier  the  abstract  factorization  theory  subsumes  the 

classical  factorization  theory  for  matrix  valued  functions  of  a complex 

variable.  It,  however,  neglects  the  question  of  structure.  For  instance, 

if  Q is  represented  by  a symetric  matrix  is  K represented  by  a symetric 

matrix  or  if  Q is  represented  by  a rational  matrix  is  K represented  by  a 

rational  matrix.  Although  the  answer  to  the  latter  question  is  known  to 

1 5 

be  yes  by  classical  factorization  theory  we  have  recently  shown  that  if 

Q is  represented  by  a rational  matrix  in  two  complex  variables  the  answer 
la 

is  no.  Indeed,  for  rational  functions  the  following  necessary  and 
sufficient  condition  may  be  derived  by  reformulating  well  known  results 

1 g 

on  several  complex  variables. 
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THEOREM:  Let  be  a polynomial  in  two  variab’es. 

then  f has  a factorization  in  the  form 

f(z^ , z^)  = g(z^ , Z2)h(z^ , z^) 

2 

where  g is  a polynomial  having  no  zeros  in  P ' h is  a polynomial 

Z 2 

having  no  zeros  in  Q = ((z^.z^)  eC  ; ! z-j  ill . ! ^2!-^  ^ 

ln| f(e^®i ,e^®2) I has  all  of  its  non-zero  Fourier  coeficients  in  the 

quadrants;  i>_0,  j^O  and  i^O,  j_<0. 

Note,  that  non-polynomial  factorizations  are  possible  under  weaker  hypotheses. 

In  essence,  the  theorem  gives  a condition  on  the  insertion  loss  function  to 

permit  it  to  be  factored  and  as  a result  a corrollary  to  the  above  theorem 

yields  a necessary  and  sufficient  condition  for  a specified  two  variable 

1 4 

insertion  loss  to  be  approximated  by  a rational  miniphase  filter. 

Relativistic  System  Theory:  Over  the  past  year  an  investigation  into  the 

possibility  of  formulating  a resolution  space  theory  in  which  a relativistic 

2 

space-time  structure  is  employed  has  been  investigated.  This  endeavor  has 
been  successful  in  the  case  of  2-dimensional  (one  space  and  one  time  dimension) 
Lorentz  space^^  though  we  have  not  been  able  to  extend  the  theory  beyond 

23 

that  case.  The  basic  difficulty  lies  with  the  need  to  generate  a semiring 

from  those  sets  in  space-time  which  are  the  futures  and  pasts  of  points. 

In  the  case  of  2-dimensional  Lorentz  space  the^ ^ "diamond  sets"  suffice  but  no 

such  semiring  apparently  exists  in  higher  dimensional  Lorentz  space  or  general 
relativistic  space-times.  In  the  case  of  2-dimensional  Lorentz  space  a resolu- 
tion space  theory  closely  paralleling  the  classical  theory  has  been  formulated 

9 

and  will  be  reported  in  a forthcoming  thesis. 

Non-Standard  Analysis:  Some  preliminary  investigations  into  the  applicability 

of  non-standard  analysis^^  to  the  discretization  of  resolution  space  have 
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been  undertaken.  Although  non-standard  analysis  is  a "logicians  game"  we 
believe  that  the  success  of  such  an  endeavor  may  lead  to  a highly  intuitive 
viewpoint  from  which  to  study  resolution  space.  In  essence,  a resolution 
space  is  a "discrete-time"  space  if  the  hermitian  operator  defined  by  the 
spectral  measure  of  the  resolution  space  has  pure  point  spectra.  When  this 
is  the  case  one  can  show  that  most  of  one's  intuition  derived  from  I2  is 
valid  in  the  abstract  space. 

Our  main  result  in  this  area  is  a new  non-standard  proof  of  a theorem 
21 

of  Berberian. 

THEOREM:  There  exists  a "faithful  functor"  from  the  catagory  of  Hil- 

bert spaces  and  bounded  linear  maps  to  itself  which  takes  each  normal 
operator  to  a normal  operator  in  which  the  entire  spectrum  is  point 
spectrum. 

By  a faithful  functor  we  mean  a functor  which  preserves  "most"  algebraic 
and  topological  Hilbert  space  properties.  Unfortunately,  the  normal  opera- 
tors of  the  theorem  do  not  have  pure  point  spectrum  and  hence  cannot  be 
used  on  our  purpose. 
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Abstract 


What  we  refer  here  to  as  "optical  noise"  actually  encompasses  a 

number  of  noise-related  problems  in  optical  information  processing. 

These  might  include  fli:i  grain  noise  and  photoelectronic  shot  noise 

in  image  processing  as  well  as  noise-induced  limitations  in  other 

applications  of  optical  information  processing.  In  these  applications 

we  must  come  to  grips  with  the  inherently  signal -dependent  nature  of 
1-2 

the  noise  sources  , in  contrast  to  the  signal-independent  noise 

3 

models  commonly  used  in  statistical  communications  problems  . It  is 

important  to  consider  the  implications  of  signal -dependent  noise  models 

in  the  detection  of  signals  in  the  presence  of  noise,  as  well  as  in 

estimating  the  parameters  of  signals.  Thus,  in  effect,  we  must  consider, 

at  the  most  basic  level,  the  implications  of  signal -dependent  noise 

models  on  the  issue  of  optimal  and  suboptimal  detectors  and  estimators 

for  applications  in  optical  information  processing.  It  is  significant 

to  note,  moreover,  that  other  (non-optical ) noise  sources,  including 

4 

magnetic  tape  recording  noise  are  effectively  signal -dependent.  Thus 
j it  appears  that  this  work  should  have  applications  to  a broad  spectrum 

of  signal -processing  problems. 

I 

Introduction 

I 

'»  To  date,  the  majority  of  work  dealing  with  s ignal -dependent  noise 

has  been  concentrated  on  rather  specialized  examples  and  applications. 

5 

I Using  a Poisson  point  process  noise  model,  Goodman  and  Belsher  have 

considered  the  restoration  of  atmospherical ly-degraded  images  using 
* linear  minimum  mean  squared  error  filters.  Walkup  and  Cheons  modified 

the  familiar  Wiener  filter  for  various  additive,  Gaussian  signal-dependent 
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2.  6 
noise  models  , and  Naderi  did  additional  work  on  this  problem®.  On 

the  other  hand.  Hunt  and  Trussel  have  derived  a nonlinear  maximum  a 

posterori  probability  (MAP)  estimate  which  can  accomodate  both  signal- 

dependent  and  signal -independent  noise  models  ' , and  have  applied  this 

MAP  estimator  to  restoring  noise-degraded  images.  For  such  applications 

and  in  the  special  case  where  the  images  of  interest  exhibit  extremely 

low  contrasts,  conventional  restoration  techniques  perform  rather  poorly. 

Therefore  hueristic  (ad  hoc)  algorithms  such  as  the  so-called  "noise 

9 

cheating"  algorithm  , have  been  developed.  Other  algorithms,  which 
explicitly  include  the  signal  dependence  of  the  noise,  as  well  as  incor- 
porating pertinent  properties  of  the  human  visual  system,  have  also  been 
investigated^.  Classical  likelihood  ratio  tests  have  also  been  derived 
for  some  rather  general  models  of  si gnal -dependent  noise^^,  but  the  solu- 
tions generally  require  a priori  knowledge  of  the  object's  exact  location, 
size,  intensity,  and  the  additive  background  present. 

In  view  of  the  rather  specific  nature  of  the  topics  investigated  to 
date,  it  is  our  opinion  that  investigations  of  a more  general,  basic 
nature  should  be  undertaken.  As  a result  we  propose  to  examine  a number 
of  rather  general  signal -dependent  noise  models  and  will  address  some 
basic  issues.  These  will  include  investigations  of  the  different  struc- 
tures of  various  optimal  estimators  designed  for  operation  in  signal- 
dependent  noise,  as  compared  with  the  structures  of  those  designed  to  work 
with  signal-indpendent  noise.  Of  particular  interest  here  are  the  issues 
of  how  and  when  a priori  information  concerning  the  "signal"  (e.g.  the 
undegraded  image  in  an  image  processing  application)  will  enable  • 
gain  additional  noise  suppression,  or  achieve  a superior  de’’-’  ’ ■ 
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Engineering  tradeoffs  relating  to  detectors  which  do,  or  do  not,  take 
the  signal-dependent  nature  of  the  noise  into  account  will  be  considered. 
Applications  to  particular  problems  in  areas  such  as  image  processing, 
optical  data  processing  with  coherent  light,  and  electronic  signal  pro- 
cessing will  be  considered  where  appropriate. 

Illustrative  Examples 

To  illustrate  the  difference  between  signal-independent  noise  pro- 
blems and  those  which  are  fundamentally  signal -dependent  in  nature,  we 
will  present  some  examples  of  particular  optimal  estimators  for  rather 
simplified  noise  models.  Consider,  for  example  the  following  two  cases: 

Case  (I)  r = s + n, 
and 

Case  (II)  r = s + ks^n, 

where  r represents  the  received  signal,  s represents  the  transmitted 

signal,  and  n represents  a noise  process.  Case  (I)  is  a familiar  textbook 

example  in  statistical  communication  theory  courses  . Case  (II)  comes 

from  image  processing,  where  the  observable  quantity  is  photographic 
2 

density  . The  constant  k might,  for  example,  be  a scanning  constant  which 
would  be  related  to  the  ratio  of  the  area  of  an  image  scanner's  aperture 
to  that  of  an  "average"  film  grain.  In  this  example  we'll  assume  that  s 
and  n are  statistically  independent,  and  that  n is  Gaussian  distributed 
with  zero  mean  and  unit  variance.  Since  in  Case  (II)  the  "noise"  is  the 
term  ks^n,  it  is  clear  that  the  noise  is  signal -dependent. 

The  effect  of  the  signal  dependence  of  the  noise  in  Case  (II)  is 
readily  apparent.  The  conditional  distribution  of  r given  s is  Gaussian 
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with  mean  s and  unit  variance  in  Case  (I),  whereas  Case  (II)  yields  a 

2 

conditional  distribution  which  is  normal  with  mean  s and  variance  k s. 

Let's  consider  the  effect  of  the  signal -dependent  nature  of  the 
noise  in  Case  (II)  on  the  structure  of  the  classical  maximum  likelihood 
estimator  (MLE).  Maximizing  the  conditional  pdf  p(r/s)  over  s,  we  find 
in  Case  (I)  that 

A 

s = r (1 ) 

where  the  circumflex  denotes  our  estimate.  Maximizing  p(r/s)  for  Case 
(II),  however,  we  find  that 

A ^ ^2  2 h .2 

s = [r^  + (|  ) ] - . (2) 

where  a negative  root  was  discarded. 

This  estimator  is  strikingly  different  from  that  of  Eq.  (1).  Note, 

however,  that  as  k-»0,  Eq.  (2)  approaches  Eq.  (1). 

Now  consider  a random  signal;  for  example  let  s be  distributed  normally 

2 

with  mean  u and  variance  o . The  appropriate  optimal  estimator  here  is  the 
maximum  a posteriori  probability  (MAP)  estimator.  Maximizing  p(s/r)  we  find 
in  Case  (I)  that 

A 2 

s = ^ ^ — r.  (3) 

a^+1  a^+1 

2 

Note  that  as  s becomes  fixed  (i.e.,  o -►0),  our  MAP  estimate  becomes 

A 

s = p , (4) 

which  ignores  entirely  the  received  signal  r.  On  the  other  hand,  as  p(s) 

2 

becomes  increasingly  disperse  (i.e.,  a -«»),  we  find  that  Eq.  (3)  reduces 
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to  Eq.  (1),  the  signal  independent  MLE.  Maximizing  p(s/r)  for  Case 
(II),  we  find  our  MAP  estimate  to  be  a solution  of  the  cubic  equation 


s'. 


- u 


2k" 


aV 

2k^ 


0. 


(5). 


Again  note  that  as  s becomes  fixed  the  solution  of  Eq.  (5)  is 

given  by  Eq.  (4)  which  again  ignores  the  received  signal  r.  Also,  as 
p(s)  becomes  disperse  (a  -+«>),  the  solution  of  Eq.  (5)  is  given  by  Eq.  (2), 
the  signal  dependent  MLE. 

As  another  example  of  the  effect  of  signal  dependent  noise,  consider 
the  familiar  Cramer-Rao  bounds.  If  ^ is  any  unbiased  estimate  of  s,  then^ 


Var  [^-s]  > {-E  - 1 


(6) 


3S 


Evaluating  Eq.  (6)  for  the  two  cases  above,  we  find  that  for  Case  (I) 

Var  [^-s]  > 1 , 


(7) 


whereas  for  Case  (II)  we  have  the  result 


Var  [^-s]  > 

2s  +k 


2 • 


(8) 


2 2s 

For  s >>  k , the  quantity  is  large  compared  to  unity,  and  the  right- 
hand  side  of  Eq.  (8)  can  be  approximated  by 


2k's'  k^s  . 


(9) 


2s+k" 


2 


In  image  processing  applications,  k is  often  much  smaller  than  unity,  so 


the  condition  s » ^ is  not  very  restrictive  ; thus,  except  for  the  very 
small  signal  case,  we  have 
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(10) 


Var  [^-s]  ^ k^s  . 

Note  that  the  lower  bound  given  by  Eq.  (10)  can  be  smaller  than  unity, 
the  lower  bound  of  Eq.  (7).  In  other  words,  the  minimum  variance  of 
our  unbiased  estimate  of  s might  actually  be  smaller  in  the  signal- 
dependent  noise  case  than  in  the  signal -independent  noise  situation. 
Preliminary  results  such  as  these  provide  ample  motivation  for  further 
investigations  into  these  problem  areas. 

Progress  on  Image  "Clutter"/"Contrast"  Investigation 

One  other  subproject  of  interest  in  the  optical  noise  area  is  the 
investigation  of  quantitative  measures  of  the  notions  of  "contrast"  and 
"clutter"  as  applied  to  imagery.  During  the  grant  period  we  have  been 
developing  the  hardware  and  software  capability  (based  on  usage  of  an 
image  storage  tube  interfaced  to  a minicomputer)  of  getting  imagery  into 
the  computer.  Our  long-range  goal  is  to  test  various  mathematical 
measures  of  contrast  and  clutter  on  the  images,  and  determine,  where 
possible,  those  measures  which  are  highly  correlated  with  the  subjective 
notions  of  image  clutter  and  contrast.  Since  these  notions  are  related  to 
a definite  form  of  psychoviaial  "noise"  in  imagery,  it  is  an  appropriate 
subject  for  study  under  the  heading  of  "optical  noise".  Presumably  one 
application  of  such  measures  of  contrast  and  clutter  would  be  in  the  area 
of  evaluating  target  location  and  identification  assistance  techniques. 

At  present  the  operating  software  is  85%  complete  and  the  system  is 
awaiting  calibration  before  the  entire  system  is  ready  for  use. 
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Abstract : 


In  the  following,  we  present  a general  model  in  which  many  questions 
arising  in  pattern  recognition  may  conveniently  be  phrased.  The  intention 
here  is  to  arrive  at  a realistic  model  which  will  allow  very  powerful 
mathematical  machinery  to  be  utilized.  The  central  thrust,  however,  of 
this  work  is  to  develop  a general  theory  of  invariance  sufficient  to  per- 
mit recognition  over  a broad  range  of  deformations  from  a standard  proto- 
type of  the  target  pattern. 

Introduction 

The  following  symbolism  will  be  used  consistently: 
n - the  set  of  patterns; 

G - a group  of  transformations  acting  on  the  left  of  n; 

K - a field  of  scalars,  usually  the  real  or  complex 
numbers; 

V - a finite  dimensional  vector  space  over  K; 

R - a map  R:  n -*■  V,  the  measurement  vector. 

In  general,  if  X and  Y are  sets  and  G acts  on  the  left  of  X,  then 
for  each  f:  X -►  Y and  each  g e G we  may  define  a new  function  gf:  X Y 
by  the  formula 

(gf)  (x)  = f(g"^x)  , xeX  . 

In  this  fashion  G acts  on  the  left  of  the  set  (F(X,Y))  of  functions 
from  X to  Y.  Note  that  G need  not  act  on  Y for  this  construction. 

With  the  above  facts  and  notation  we  may  present  the  basic  model. 

The  basic  idea  is  that  the  natural  transformations  of  the  patterns  are 
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assumed  to  be  well  modelled  by  the  action  of  elements  of  the  trans- 
formation group  G.  Thus,  if  w e n and  g e G,  then  gw  is  another 
pattern  which  has  been  obtained  as  a transform  of  w.  Finally,  we  per- 
form certain  measurements  on  w,  each  of  which  results  in  a scalar  value 
X e K.  The  collection  of  all  such  measurements  is  assembled  as  a 
measurement  vector,  R(w)  e V. 

Now,  we  may  give  two  interpretations.  For  the  first  of  these,  we 
regard  R as  a fixed  retina  or  camera  and  think  of  the  group  of  trans- 
formations as  moving  the  patterns  under  the  camera.  In  the  second 
interpretation,  through  the  induced  action  of  G on  F(fi,V)  discussed 
above,  we  imagine  the  group  as  acting  on  the  camera.  Thus,  the  group 
G can  serve  as  a parameter  set  for  the  control  mechanism  associated 
with  the  measurement  apparatus.  It  is  important  that  we  observe  that 
the  mathematical  formalism  is  exactly  the  same  irrespective  of  which 
interpretation  is  applicable  to  a particular  problem.  In  fact,  many 
situations  involve  some  combination  of  both  concepts. 


In  summary,  we  have  represented  as  a set  of  functions,  each  from 

G to  V.  Accordingly,  we  will  refer  to  a map  R : n V as  a representa- 
tion of  Q in  V. 

We  have  constructed  a map  w (->  w from  n into  F(G,  V),  Let  us  note 
that  the  following  property  holds: 

g w = fg^,  g £ G,  w e n, 

where  g w is  obtained  from  the  action  of  G on  F(G,  V).  This  leads  us  to 
introduce  the  following: 

Definition:  A functional  representation  of  Q in  a class  F of  functions 

from  G to  V is  a map  r : n -f  F,  denoted  w i-<-  w"^,  such  that 

g w*"  = (gw)'" 

for  all  g E G and  w e a. 

The  representation  w h-  w induced  by  R : n ->■  V is  an  example  of  a 
functional  representation.  We  can  show  that  this  is  not  coincidental. 

In  fact,  we  have  the  following: 

Theorem:  The  functional  representations  r of  n as  functions  frcm  G to  V 

correspond  one-to-one  with  the  representations  R:  n-*V  via  the  corres- 
pondence  r ♦-*■  R if  and  only  if  w (g)  = R(g  w)  for  all  g e G and 
wen. 

The  connection  established  by  the  preceding  theorem  is  a very  natural 
one  and  has  significant  applications.  The  maps  R and  r of  the  Theorem 
will  henceforth  be  assumed  related  as  above. 

In  effect,  this  result  permits  the  transfer  of  structure  from  the 
group  G (more  or  less  concrete)  to  the  set  n of  patterns  (an  abstract 
entity).  For  example,  let  us  suppose  that  G is  a topological  group  ^ 

I 
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and  that  K is  the  real  or  comples  numbers . We  define  a representation 
R : r,  -*■  V to  be  continuous  provided  that  each  w c fi,  is  continuous 
in  the  usual  sense.  Nov;  let  Cy(o)  denote  the  set  of  representations 
which  are  continuous  in  this  sense.  It  is  reasonable  to  ask  if  there 
exists  a topology  T on  the  set  such  that  C^(n)  is  precisely  the  class 
of  T-  continuous  functions  from  G to  V,  i.e.,  that  C^{2)  = C^(q,T). 

This  is  indeed  the  case,  and  the  weakest  such  T has  useful  properties 
with  respect  to  continuity  of  the  action  of  G on  It  is  quite  likely 
that  our  choice  of  a camera  should  be  restricted  to  continuous  repre- 
sentations. This  merits  additional  consideration. 

9 

In  general,  the  representation  theorem  above,  together  with  the 
overall  smoothing  effect  normally  found  in  physical  measurements  justify 
making  mathematically  useful  assumptions  about  the  set  of  patterns.  For 
instance,  we  might  begin  a particular  investigation  with  "By  a pattern 
we  mean  a continuous  function  from  a topological  group  G into  a finite 
dimensional  vector  space  V ...".  Likewise,  continuity  can  be  replaced 
by  a variety  of  other  properties  such  as  Haar  integrable,  differentiable, 
etc. 

Formulation  of  Two  Classical  Problems 

I . Template  Matching 

The  idea  here  is  to  compare  a pattern  against  a number  (usually 
finite)  of  prototypes  for  goodness  of  fit.  The  problem  is  well  worked 
over  in  the  literature  in  many  situations.  However,  the  work  is  some- 
what limited  in  the  case  in  which  t’.e  object  is  to  match  the  prototypes 
against  a possible  transform  of  the  pattern.  Many  investigations  do 
incorporate  a limited  amount  of  correction  for  such  effects  as  translation 
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and  magnifications.  However,  rotational  effects  seem  to  be  generally 
avoided. 

In  the  context  of  our  model,  where  a choice  of  representation  R 
allows  us  to  deal  with  the  corresponding  functional  representation, 
the  question  of  matching  becomes;  given  w^ , W2  e a,  find  g e G such 
that  g w^  = w^,  i.e.,  w^(g  ^x)  = W2(x)  for  all  x e G.  Then  questions 
of  goodness  of  fit  become  questions  of  approximation  theory. 

We  anticipate  results  of  two  types  on  the  template  matching  pro- 
blem. First,  we  should  be  able  to  determine  the  existence  or  non- 
existence of  a solution  for  g e G in  the  equation  gw^  = W2*  Hopefully, 
we  will  be  able  to  develop  invariant  representations  or  relatively  in- 
variant representations  for  which  the  solution  to  the  matching  problem 
is  straightforward.  Finally,  we  whould  be  able  to  show  that,  in  a 
random  environment,  the  use  of  Gaussian  statistics  can  be  replaced  by 

purely  analytic  methods  involving  least  squares  approximation  in  an 
2 

L space. 

The  major  tool  here  will  be  the  Haar  integral  and  the  assumptions 
necessary  for  pursuit  along  the  above  lines  seem  to  be  rather  weak. 
Perhaps  a locally  compact  group  is  sufficient. 

II . Imbedded  Patterns 

The  idea  here  is  to  search  a pattern  for  the  local  presence 
of  a predetermined  sub-pattern. 

This  problem  appears  to  be  one  which  gave  much  impetus  to  pat- 
tern recognition  and  artifical  intelligence  during  the  early  1960's. 
However,  anticipated  success  has  failed  to  materialize.  Perhaps  this  is 
due  to  the  failure  to  introduce  enough  powerful  mathematical  techniques. 
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In  the  present  context  we  may  formulate  a simple  version  of  this 
problem  as  follows;  Given  w c a and  v^  c V,  find  g e G (if  possible) 
such  that  w (g)  = V^.  Thus,  the  problem  reduces  to  the  solution  of 
a functional  equation.  An  analytic  structure  of  some  type  is  in  order 
for  this  problem,  and  we  propose  the  use  of  Lie  groups.  Upon  intro- 
duction of  coordinates  in  G we  obtain  a system  of  non-linear  equations 
and  seek  a "best"  solution. 

An  attempt  to  minimize  the  objective  function 

f (x)  = II  w'"(x)  - 11^/2 

yields  the  differential  equation 

X (t)  = A*(x)  [w'"(x)-Vq]  , x(0)=Xq, 

where  x e k'*’  is  the  coordinate  vector  of  a group  element  and  A(x)  is 
the  derivative  of  w (x).  This  equation  can  be  solved  by  ordinary 
numerical  methods,  although  we  need  to  develop  sophisticated  methods 
which  utilize  the  group  structure  and  avoid  the  differentiation  of 
w^(x).  We  fine  that  x = lim  x(t)  gives  a critical  value  for  the  ob- 

t-H» 

jective  function,  'i'(x). 

It  is  hopeful  that  some  combination  of  techniques  such  as  hill- 
climbing together  with  the  solution  of  differential  equations  will  pro- 
vide a fast  and  efficient  approach  to  the  solution  of  this  problem. 

Summary 

We  hope  to  further  develop  the  general  framework  for  certain  pattern 
recognition  problems  involving  transformation  groups.  Additional  work 
is  to  be  done  with  regard  to  the  transfer  to  topological  and  analytical 
structure  from  the  group  to  the  pattern  space.  Through  the  use  of  in- 
variants, relative  invariants,  Haar  integration,  and  the  solution  of 
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problems  of  extrema  we  hope  to  obtain  workable  results  for  use  in  con- 
nection with  the  template  matching  problem  and  the  imbedded  pattern 
problem. 

The  immediate  interest  will  involve  compact  groups.  This  is 
partly  because  the  representation  theory  for  such  groups  is  known  in 
good  detail  and  also  because  it  includes  an  important  special  case, 
namely  the  rotation  groups,  which  has  not  received  an  adequate  treat- 
ment in  the  literature. 
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742-3506 

R.  Saeks 

Assoc.  Prof. -EE 
and  Math 

Circuits  and 
Systems 

258A-EE 

742-3528 

J.  Walkup 

Assoc.  Prof. -EE 

Optical  Systems  & 
Communications 

260B-EE 

742-3500 

D.  Vines 

Prof. -EE 

Computers 

151-EE 

742-3536 

PHYSICAL  ELECTRONICS 

K.  DasGupta 

Prof .-Phys. 

X-Ray  Spectroscopy 

4-Sci . 

742-3773 

M.  Gundersen 

Asst.  Prof. -EE 

Quantum  Electronics 

260A-EE 

742-3501 

W.  Portnoy 

Prof. -EE 

Sol  id  State  & Bio- 
Medical  Electronics 

152-EE 

742-3532 

R.  Redington 

Prof .-Chem. 

Molecular  Spectro. 

27-Chem 

742-3088 

J.  Reichert 

Assoc.  Prof-EE 

Optics  & Quantum  Elec 

203-EE 

742-3502 

G.  Robinson 

Prof . -Chem. 

Transient  Phenomena 

313A-Chem 

742-3095 

H.  Thomas 

Prof .-Phys. 

Solid  State 

75-Sci. 

742-3780 

R,  Wilde 

Asso.  Prof-Chem 

Molecular  Spectro. 

26-Chem 

742-3094 

F.  Williams 

Asst.  Prof-EE 

Interaction  of  Light 
with  Matter 

258B-EE 

742-3501 
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ELECTROMAGNETICS 


Name 

Title 

Area 

Office 

Phone 

J.  Craig 

Prof. -EE 

Radiation  Analysis 

101-EE 

742-3529 

R.  Cross 

Vis.  Prof-EE 

Plasma 

103C-EE 

742-3468 

M.  Hagler 

Prof. -EE 

Plasma 

103B-EE 

742-3470 

M.  Kristiansen 

Prof. -EE 

Plasma 

103A-EE 

742-3468 

E.  Kunhardt 

Asst.  Prof-EE 

Nonlinear  Phenomena 

260C-EE 

742-3545 

T.  Trost 

Assoc.  Prof-EE 

Antennas 

102-EE 

742-3505 

POWER 

T.  Burkes 

Assoc.  Prof-EE 

Power  Conditioning 

105C-EE 

742-3533 

J.  Craig 

Prof-EE 

Electro-mech.  devices 

101-EE 

742-3529 

E.  Kunhardt 

Asst.  Prof-EE 

High  Power  Switching 

260C-EE 

742-3545 

S.  Liberty 

Assoc.  Prof-EE 

Solar  Energy 

201B-EE 

742-3441 

J.  Reichert 

Assoc.  Prof-EE 

Solar  Energy 

203-EE 

742-3502 

T.  Trost 

Assoc.  Prof-EE 

Magnetic  Energy  Strg 

102-EE 

742-3505 

1 

¥ 


i 


» 


ACTIVE  GRANTS  AND  CONTRACTS 
in 

ELECTRONICS  AND  RELATED  AREAS 


•I 

I 


) 


* - J 


Systems 


» 


* 


9 


Principal 

Invest. 

Agency 

Title 

Duration 

Total 

Funding 

Present 
Annual  Fund 

Saeks 

AFOSR 

Resolution  Space, . . 

12/73-4/78 

$ 80,800 

$ 23,500 

Liberty 

ONR 

Kalman  Filtering 

3/74-3/78 

93,300 

25,000 

Saeks 

ONR 

Fault  Analysis 

4/74-12/77 

113,000 

30,000 

Walkup/ 

Hagler 

NSF 

Optical  Info.  Proc. 

6/75-1/77 

20,600 

14,000 

Chao/Saeks 

NSF 

Analysis  & Design. . 

6/75-11/77 

30,000 

15,000 

Wal kup/ 
Hagler 

AFOSR 

Volume  Hologram/ 
Optical  Syst. 

6/74-9/77 

133,000 

67,000 

Saeks 

ONR 

Unified  Prg.  in  Elec. 

.9/76-2/78 

200,000 

133,000 

Gustafson 

SORF 

Microprocessors 

9/76-8/77 

4,500 

4,500 

Total  Annual 
in  Systems 

Funding 

$312,500 

Physical  Electronics 


■B 

Title 

■IQnm 

Present 
Annual  Fund 

Gundersen 

ERDA 

Pulsed  Molecular. . 

7/74-9/77 

$191,000 

$ 75,000 

Gundersen 

NSF 

Infrared  Upconversion 

11/75-10/77 

35,000 

18,000 

Reichert 

NSF 

Undergraduate  Research  3/71  2/77 

90,000 

16,000 

Portnoy 

NASA 

Insulated  ECG  Elect. . 

7/71-277 

85,000 

25,000 

Reichert 

AFOSR 

Unstable  Opt.  Res. 

10/72-12/76 

200,000 

50,000 

Portnoy 

SORF 

Solid  State  Research 

9/76-8/77 

5,000 

5,000 

Williams 

SORF 

Quantum  Electronics 

9/76-8/77 

5,000 

5,000 

Thomas 

Welch 

Low  Temp  Solid  State 

6/75-5/78 

45,000 

15,000 

Robinson 

ARO 

Raman  Scattering 

9/76-8/77 

25,000 

25,000 

Robinson 

Welch 

Transient  Studies 

9/76-8/79 

150,000 

50,000 

Wilde 

Welch 

Sol  id  State  Studies 

6/76-5/79 

51 ,000 

17,000 

Total  Annual  Funding 
in  Physical  Electronics  $ 301,000 


A 
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Electromagnetics 


Principal 

Invest. 

Agency 

Title 

Duration 

Total 

Funding 

Present 
Annual  Fund 

Kristiansen 

AFOSR 

Plasma  Heat  for  Therm 

3/71-10/77 

$410,000 

$100,000 

Kristiansen 

NSR 

Toroidal  Plasma  Heat. 

6/76-5/77 

24,500 

24,500 

Trost 

SNF 

Radio  Bursts,  Tornado 

10/75-5/77 

42,500 

21,000 

Kristiansen 

SORF 

Toriodal  Plasma  Heat. 

6/76-5/77 

24,500 

24,500 

Kri stiansen 

NSF 

RF  Plasma  Heating 

4/71-5/77 

205,600 

30,000 

Total  Annual  Funding 
in  Electromagnetics 

$200,000 

Power 

Principal 

Invest. 

Agency 

Title 

Duration 

Total 

Funding 

Present 
Annual  Fund 

Craig 

TPL 

Power  Systems 

1/73-12/76 

$ 32,000 

$ 8,000 

Burkes 

AFWL 

Airborne  Power  Systs. 

2/74-5/77 

167,300 

57,000 

Trost 

NSWC 

Magnetic  Energy  Strg. 

12/75-12/76 

16,645 

16,645 

Burkes 

ERDA 

E Beam  Laser 

3/76-1/77 

23,000 

23,000 

Burkes 

AFOSR 

Pulsed  Power  Conf. 

6/76-5/77 

6,000 

6,000 

Kristinasen 

AFWL 

High  Power  Switch  Dev 

9/76-9/77 

50,000 

50,000 

Reichert/ 

ERDA 

Crosbyton  Solar  Power 

9/76-7/77 

880,000 

400,000 

Liberty 

♦Includes  sub-contractors  funding 


Total  Annual  Funding 
in  Power  $560,645 


Total  Annual  Funding 
in  Electronics  and 
Related  Areas  $1 ,374,145 


RESEARCH  LABORATORIES 
in 

ELECTRONICS  AND  RELATED  AREAS 


SYSTEMS 


Computer  Laboratories: 

CDC  1604  facility:  hands-on  facility  for  both  education 

and  research 108-EE 

Hybrid  Computer  facility:  minis,  micros  and  analog 

facil  ities 162-EE 

Bio-medical  Systems:  includes  instrumentation  and  microprocessor 

application  facilities 208-EE 

Circuits  and  Systems  Laboratory:  the  think  tank 258-EE 

Optical  Systems  Laboratories: 

Holographic  Optics:  primarily  used  for  multiplex  holography 

research 11 0-EE 

Optical  Signal  Processing:  research  in  optical  and  digital 

image  processing 216-EE 

PHYSICAL  ELECTRONICS 

Laser  Laboratory:  infrared  laser  research ...262-EE 

Solid  State  Laboratories: 

Solid  State  Physics:  basic  semiconductor  studies 34-Sci. 

Low  Temperature  Physics:  superconductor  studies 58-Sci. 

Integrated  Circuit:  fabrication  facility  for  SSI  and 

special  purpose  devices 209-EE 

Spectroscopy  Laboratories: 

Laser  Spectroscopy:  interaction  of  light  with  matter 260-EE 

Molecular  Spectroscopy:  semiconductor  studies 35-Chem 

Picosecond  Spectroscopy:  transient  studies 30-Chem 

X-Ray  Spectroscopy:  X-ray  studies  and  laser  development 2-Sci. 

ELECTROMAGNETICS 
Plasma  Laboratories: 

Laser/Plasma  facility:  plasma  heating  via  laser  plasma 

interaction 113-EE 
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ELECTROMAGNETICS  (continued) 


Plasma  Laboratories  (continued); 

Tokamak  facility:  radio  frequency  heating  of  toroidal 


plasmas 1 17-EE 

Electromagnetics  Laboratory:  nonlinear  wave  studies 111-EE 


Antenna  Laboratory:  radio  meteorology  and  ionspheric 

studies West  of  the 

Medical  School 


Power 

High  Voltage  Laboratory;  pulsed  power  studies. . .North  of  Textile  Bldg 


Solar  Energy  Laboratory:  another  think  tank 205-EE 

High  Power  Switching  Laboratory:  electron  beam  initiated 

spark  gap Trailer  west  of  EE 

Bldg. 


PUBLICATION  ACTIVITY 
in 

ELECTRONICS  AND  RELATED  AREAS 


Systems  - Refereed  Publications 


Liberty,  S.,  "On  the  Essential  Quadratic  Nature  of  LQG  Control  - 
Performance  Measure  Cumulants,"  Information  and  Control,  Vol . 32, 

No.  3,  pp.  276-305,  Nov.  1976. 

Saeks,  R.,  Large  Scale  Dynamical  Systems  No.  Hollywood,  Point  Lobos 
Press,  1976^  (314  + v pages,  edited  volume  of  contributed  chapters). 

Saeks,  R.,  "Analysis  and  Design  of  Interconnected  Dynamical  Systems," 
in  Large-Scale  Dynamical  Systems,  No.  Hollywood,  Point  Lobos  Press, 
pp.  59-79,  (with  G.  Wise  and  K.S.  Chao). 

Saeks,  R.,  "The  Factorization  Problem  - A Survey,"  IEEE  Proc.,  Vol. 

64,  pp.  90-95,  (Jan.  1976). 

Saeks,  R.,  "Continuations  Approach  to  Large-Change  Sensitivity  Analysis," 
IEEE  (London)  Journal  on  Electronic  Circuits  and  Systems,  Vol.  1,  pp. 

11-15  I'Oct.  1976,  with  K.S.  Chao). 

Saeks,  R. , "A  Spectral  Theoretic  Approach  to  Fault  Analysis  in  Linear 
Sequential  Circuits,"  Journal  of  the  Franklin  Inst.,  Vol.  302,  pp. 

239-258,  (with  S.  R.  Liberty  and  S.  Sangani , Sept.  1976). 

Saeks,  R.,  "Representation  of  Weakly  Additive  Operators,"  Proc.  of  the 
A.M.S.,  Vol.  59,  pp.  55-61,  (with  R.A.  DeCarlo,  Aug.  1976). 

Hagler,  M. , "A  Sampling  Theorem  for  Space-Variant  Systems,”  J.  Opt.  Soc. 
Am.  66,  918  (1976),  (with  R.  J.  Marks  II  and  J.F.  Walkup). 

Hagler,  M.,  "Line  Spread  Function  Notation,"  Appl . Optics  1_5,  2289  (1976), 
(with  R.  J.  Marks  II  and  J.F.  Walkup). 

Systems  - Conference  Papers  and  Reports 

Gustafson,  D.,  "A  Microprogrammed  Machine  Architecture  for  Efficient 
Matrix  Multiplication,"  SIGMICRO  Vol.  7 No.  3,  Sept.  1976,  pp.  56-61. 

Walkup,  J.,  "An  Improved  Coherent  Processor  for  Ambiguity  Function  Dis- 
play," Int.  Optical  Computing  Conf . , Capri,  Italy,  Sept.  1976,  (Proc.  in 
press, IEEE  Press)  (with  R.  J.  Marks  II  and  T.  F.  Krile). 

Liberty,  S.  R. , "A  New  Class  of  Statistical  Performance  Criteria  for 
Stochastic  Linear  Control  Systems,  " Proc.  IEEE  Conf.  on  Decision  and 
Control,  pp.  1136-1143,  Dec.  1976. 

Portnoy,  W.  M.,  "A  Recognition  System  for  the  Detection  of  Cardiac  Arr- 
Tiytanias,"  Proc.  of  the  IEEE  Conference  on  the  Applications  of  Electronics 
in  Medicine,  Number  34,  223  (1976). 

Prabhakar,  J.,  "Digital  Filters  Using  Walsh  Functions,"  IEEE  Region  V 
Conf.  April  1976,  Austin,  Texas. 

Vines,  D.,  "Automation:  An  Electrical  Engineers  Perspective,"  Rocky  Moun- 

tain Regional  ASCS  Conference,  Wyoming. 
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Systems  - Conference  Papers  and  Reports  (continued) 

Saeks,  R. , "Fault  Analysis  in  Affine  Sequential  Circuits,"  Proc.  of 
the  1976  Conf.  on  Info.  Sci.  and  Systems,  pp.  227-232,  Johns  Hopkins 
Univ.,  March  1976  (with  S.  Sangani). 

Saeks,  R. , "On  the  Computation  of  Rational  z-Transforms ,"  Proc.  of  the 
1976  IEEE  Inter.  Symp.  on  Circuits  and  Systems,  pp.  518-521,  Munich, 

April  1976  (with  S.  Sangani  and  S.  R.  Liberty). 

Saeks,  R.,  "An  Experimentation  in  Fault  Prediction  II,"  Proc.  of  the 
IEEE  AUTOTESTCON,  Nov.  1977,  Arlington,  Texas,  p.  53  (abstract  only). 

Saeks,  R. , "A  New  Characterization  of  the  Nyquist  Stability  Criterion," 
Proc.  of  the  19th  Midwest  Symp.  on  Circuits  and  Systems,  Univ.  of  Wis- 
consin at  Milwaukee,  pp.  349-353,  Aug.  1976,  (with  R.  A.  DeCarlo). 

Chao,  K.S.,  "Continuation  Methods  for  Stability  Analysis  of  Multi- 
variable  Feedback  Systems,  "Proc.  19th  Midwest  Symposium  on  Circuits 
and  Systems  (with  E.  Huang  and  R.  Saeks). 

Hagler,  M.,  "General  One-Dimensional  Space-Variant  Coherent  Processors," 
1976  Meeting,  Optical  Soc.  of  America,  October  1976  (with  R.  J.  Marks  II 
and  J . F.  Wal kup) . 

Hagler,  M.,  "Optical  Inform.ation  Processing  Experiments  for  Undergraduate 
Engineers,"  1976  Meeting,  Optical  Soc.  of  America,  October  1976  (with  G.  R 
Forehlich  and  J.F.  Walkup). 

Hagler,  M.,  "A  Spatial  Signal  Feedback  Amplifier,"  1976  Meeting,  Optical 
Society  of  America,  October  1976  (with  S.  V.  Bell). 

Systems  - in  press 

Walkup,  J.F.,  "Ambiguity  Function  Display:  An  Improved  Coherent  Processor 

Appl . Optics.  , in  press,  (with  R.  J.  Marks  II  and  T.  F.  Krile). 

Portnoy,  W.,  "A  Nonlinear  Black-Box  Electrical  Model  of  a Mesquite  Wood 
Thermal  System,"  Journal  of  the  Franklin  Institute  (in  press). 

Saeks,  R. , Interconnected  Dynamical  Systems,  New  York,  Marcel  Dekker,  Inc. 
(with  R.  A.  DeCarlo,  in  press). 

Saeks,  R.,  "Reproducing  Kernel  Resolution  Space  and  its  Applications," 
Jour,  of  the  Franklin  Institute  (to  appear). 

Saeks,  R. , "Fault  Prediction  - Twoards  a Mathematical  Theory,"  in  Rational 
Fault  Analysis,  New  York,  Marcel  Dekker,  Inc.  (with  L.  Tung  and  S.  R. 
Liberty,  to  appear). 

Saeks,  R. , "A  Functional  Approach  to  Fault  Analysis,"  in  Rational  Fault 
Analysis,  New  York,  Marcel  Dekker,  Inc.  (with  M.N.  Ramsom,  to  appear) . 
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Systems  - in  press  (continued) 


Seeks,  R.  , "An  Algebraic  Topological  Proof  of  the  Nyquist  Criterion," 

Inter.  Jour,  on  Control  (with  R.  A.  DeCarlo,  to  appear). 

Seeks,  R. , "Multivariable  Nyquist  Theory,"  Inter.  Jour,  on  Control, 

(with  R.  A.  DeCarlo  and  J.  Murray,  to  appear). 

Seeks,  R. , "A  Nyquist-like  Stability  Test  for  Two  Variable  Digital 
Filters,"  IEEE  Proc.,  (with  R.  A.  DeCarlo  and  J.  Murray,  to  appear). 

Seeks,  R. , "Three  Graphical  Stability  Tests  for  Two-Dimensional  Digital 
Filters,"  Proc.  of  the  1977  IEEE  Inter.  Symp.  on  Circuits  and  Systems, 
Phoenix,  April  1977  (with  R.  A.  DeCarlo  and  J.  Murray,  to  appear). 

Seeks,  R. , "On  the  Computability  of  Lacunary  Sets,"  Proc.  of  the  1977 
IEEE  Inter.  Symp.  on  Circuits  and  Systems,  Phoenix,  April,  1977  (with 
R.  J.  Leake  and  W.  J.  Kelly,  to  appear). 

Saeks,  R.  , "Fault  Analysis  via  Affinization,"  Proc.  of  the  1976  Asilomar 
Conf.  on  Circuits,  Systems  and  Computers,  Pacific  Grove,  CA. , Nov.  1976 
(with  S.  Sangani  and  R.  A.  DeCarlo,  to  appear). 

Chao,  K.S.,  "Transfer-Characteristic  Plots  and  Large-Change  Sensitivity 
of  Nonlinear  Resistive  Networks,"  Symposium  on  Circuits  and  Systems, 

(with  D.  K.  Liu,  in  press). 

Chao,  K.S.,  "Continuation  Methods  in  Circuit  Analysis,"  Proc.  of  IEEE, 

(with  R.  Saeks,  in  press). 

Hagler,  M. , "Space-Variant  Processing  of  One-Dimensional  Signals,"  Appl . 
Optics,  (with  R.  J.  Marks  II,  J.  F.  Walkup  and  T.  F.  Kril,  in  press). 

Hagler,  M.,  "Volume  Hologram  Representation  of  Space-Variant  Systems," 
Inter.  Conf.  on  Applications  of  Holography  and  Optical  Data  Processing, 
Jerusalem,  Israel,  August,  l''-76  (proc.  in  press,  Pergamon  Press  Ltd.) 

(with  R.  J.  Marks  II  and  J.  F.  Walkup). 

Hagler,  M. , "The  Design  of  a Spatial  Signal  Feedback  Amplifier,"  Int. 
Optical  Computing  Conference,  Copier,  Italy,  Sept.  1976  (proc.  in  press, 
IEEE  Press)  (with  S.  V.  Bell). 

Physical  Electronics  - Refereed  Publications 

Gundersen,  M.,  "Effect  of  a Small  Capacitor  in  Parallel  with  a Pulsed  COo 
TEA  Laser,"  IEEE  J.  Quantum  Electron.  QE-1 2 , 447  ( 1976)  (with  A.  H.  Bush- 
nel 1 and  T.  R.  Burkes) . 

1 5 

Gundersen,  M.,  "New  Far-Infrared  Lines  in  N H,,"  IEEE  J.  Quantum  Electron. 
QE-12,  260  (1976)  (with  A.  H.  Bushnell). 

Portnoy,  W.,  "Technology  in  Art:  An  Undergraduate  Experience,"  Inter. 

Journal  of  Electrical  Engineering  Education  ]_3,  206  (1976). 
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Physical  Electronics  - Refereed  Publication  (continued) 

Williams,  P.  F.,  "Some  Comments  on  the  Plasmon  Spectrum  of  Tetrathiaful - 
valene  Tetracyano-p-quinodimethane  (TTF-TCNQ),"  Phys,  Rev.  Lett  36,  64 
(1976)  (with  A.  N.  Bloch). 

Williams,  P.F.,  "Resonance  Raman  Scattering  of  Light  from  a Diatomic 
Molecule,"  J.  Chem.  Phys.  3519  (1976)  (with  D.  L.  Rousseau). 

Williams,  P.F.,  "Comment  on  Two-Magnon  Resonant  Raman  Scattering  in  MnFp," 
Phys.  Rev.  Lett.  (with  D.  L.  Rousseau,  R.  E.  Dietz  and  H.  J.  Guggen- 
heim) . 

Williams,  P.F.,  "Reply  to  'Comment  on'  Resonance  Raman  Scattering  and 
Collision  - Induced  Redistribution  Scattering  in  Ip,  Phys.  Rev.  Lett. 

1441  (1976)  (with  D.  L.  Rousseau  and  G.  D.  Patterson). 

Redington,  R.  L.,  "Kinetics  of  Oxygen-18  Exchange  Between  Carboxylic  Acids 
and  Water,"  J.  Phys.  Chem.,  M,  229  (1976). 

Wilde,  R.  E.,  "Detection  of  Phase  Transitions  in  Heavy  Silane  by  Infrared 
Spectroscopy,"  J.  Phys.  Chem.  Solids,  36,  119  (1975)  (with  T.K.K.  Srini- 
j vasan) . 

Wilde,  R.  E.,  "Curtius  and  Lossen  Rearrangements.  III.  Photolysis  of  Cer- 
tain Carbemoyl  Azides,"  J.  Org.  Chem.  2608  (1975)  (with  W.  Lwowski, 

R.  A.  de  Mauriac,  M.  Thompson  and  S.  Y.  Chen). 

Wilde,  R,  E.,  "An  I.R.  Study  of  Heavy  Phosphine  in  the  Solid  State,"  J. 

Phys.  Chem.  Solids,  36,  1225  (1975)  (with  B.  C.  Covington). 

Robinson,  G.  W.,  "Rhodopsin  Cooperativi ty  in  Visual  Responses,"  Vision  Res. 
15,  35  (1975). 

I Robinson,  G.  W.,  "Formation  of  Molecules  on  Small  Interstellar  Grains," 

‘ Astrophys.  J.,  195,  81  (1975)  (with  M.  Allen). 

I 

Robinson,  G.  W. , "The  Level  Shift  Operator  and  its  Effect  on  Line  Shapes 
' in  Vibronically  Perturbed  Spectra,"  Mol.  Phys.,  ^9,  613  (1975)  (with  C. 

j A.  Langhoff). 

t*  Robinson,  G.  W. , "An  Approach  to  the  Understanding  of  Radiation  Chemistry 

in  the  Condensed  Phase,"  Chem.  Phys.  Lett.,  211  (1975)  (with  J.O.  Berg). 

Robinson,  G.  W.,  "Lineshape-Lifetime  Relationship  and  Emission  and  Scatter- 
' ing  of  Light  by  Polyatomic  Molecules,"  Can.  J.  Phys.,  13,  2068  (1975) 

Herzberg  Festschrift  Edition,  (with  J.  0.  Berg). 

Thomas,  H.C.,  "Impurity  Conduction  in  Transmutation-Doped  Germanuim  During 
Decay  of  the  Radioactive  Products,"  J.  Appl.  Phys.  46  4541  (1975)  (with 
B.  Covington). 
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Physical  Electronics  - in  press 


Portnoy,  W.,  "A  Digital  Capacitance  Meter  for  Foliage  Yield  Measure- 
ments," Proc.  of  the  Seventh  Congress  of  the  International  Measurement 
Confederation  (in  press). 

Williams,  P.  F. , "Experimental  Measurement  of  Dissociative  Molecular 
Potential  Functions  from  Continuum  Resonance  Raman  Spectra,"  (with  A. 
Fernandez  and  D.  L.  Rousseau)  (accepted  for  publication  in  Chem. 

Phys.  Lett.) 

Redington,  R.  L.,  "Vibrational  Spectra  and  Normal  Coordinate  Analysis 
of  Isotopically  Labelled  Formic  Acid  Monomers,"  J.  Molec.  Spectrosc., 

( in  press) . 

Robinson,  G.  W.,  "Extraction  of  Vibronic  Information  from  Tangled 
Spectra,"  (with  J.  0.  Berg)  (Chem.  Phys.  Lett.,  in  press). 

Robinson,  G.  W.  , "Molecular  Hydrogen  in  Interstellar  Dark  Clouds, "(with 
M.  Allen)  (Astrophys.  J.,  in  press). 
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